| 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 // G4TwistTubsFlatSide implementation <<
27 // 26 //
28 // 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepbu << 27 // $Id: G4TwistTubsFlatSide.cc 72937 2013-08-14 13:20:38Z gcosmo $
29 // 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), << 28 //
30 // from original version in Jupi << 29 //
>> 30 // --------------------------------------------------------------------
>> 31 // GEANT 4 class source file
>> 32 //
>> 33 //
>> 34 // G4TwistTubsFlatSide.cc
>> 35 //
>> 36 // Author:
>> 37 // 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepburn.s.chiba-u.ac.jp)
>> 38 //
>> 39 // History:
>> 40 // 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), Integration in Geant4
>> 41 // from original version in Jupiter-2.5.02 application.
31 // ------------------------------------------- 42 // --------------------------------------------------------------------
32 43
33 #include "G4TwistTubsFlatSide.hh" 44 #include "G4TwistTubsFlatSide.hh"
34 #include "G4GeometryTolerance.hh" 45 #include "G4GeometryTolerance.hh"
35 46
36 //============================================ 47 //=====================================================================
37 //* constructors ----------------------------- 48 //* constructors ------------------------------------------------------
38 49
39 G4TwistTubsFlatSide::G4TwistTubsFlatSide(const << 50 G4TwistTubsFlatSide::G4TwistTubsFlatSide(const G4String &name,
40 const G4RotationM << 51 const G4RotationMatrix &rot,
41 const G4ThreeVect << 52 const G4ThreeVector &tlate,
42 const G4ThreeVect << 53 const G4ThreeVector &n,
43 const EAxis 54 const EAxis axis0 ,
44 const EAxis 55 const EAxis axis1 ,
45 G4double 56 G4double axis0min,
46 G4double 57 G4double axis1min,
47 G4double 58 G4double axis0max,
48 G4double 59 G4double axis1max )
49 : G4VTwistSurface(name, rot, tlate, 0, axis0 60 : G4VTwistSurface(name, rot, tlate, 0, axis0, axis1,
50 axis0min, axis1min, axis0m << 61 axis0min, axis1min, axis0max, axis1max)
51 { 62 {
52 if (axis0 == kPhi && axis1 == kRho) << 63 if (axis0 == kPhi && axis1 == kRho) {
53 { <<
54 G4Exception("G4TwistTubsFlatSide::G4Twis 64 G4Exception("G4TwistTubsFlatSide::G4TwistTubsFlatSide()",
55 "GeomSolids0002", FatalError 65 "GeomSolids0002", FatalErrorInArgument,
56 "Should swap axis0 and axis1 66 "Should swap axis0 and axis1!");
57 } 67 }
58 68
59 G4ThreeVector normal = rot.inverse()*n; 69 G4ThreeVector normal = rot.inverse()*n;
60 fCurrentNormal.normal = normal.unit(); // 70 fCurrentNormal.normal = normal.unit(); // in local coordinate system
61 fIsValidNorm = true; 71 fIsValidNorm = true;
62 72
63 SetCorners(); 73 SetCorners();
64 SetBoundaries(); 74 SetBoundaries();
65 75
66 fSurfaceArea = 1. ; // NOTE: not yet imple << 76 fSurfaceArea = 1 ; // not yet implemented. This is NOT a problem for tracking
>> 77
67 } 78 }
68 79
69 G4TwistTubsFlatSide::G4TwistTubsFlatSide( cons << 80
70 << 81
71 << 82 G4TwistTubsFlatSide::G4TwistTubsFlatSide( const G4String &name,
72 << 83 G4double EndInnerRadius[2],
73 << 84 G4double EndOuterRadius[2],
74 << 85 G4double DPhi,
75 << 86 G4double EndPhi[2],
>> 87 G4double EndZ[2],
>> 88 G4int handedness )
76 : G4VTwistSurface(name) 89 : G4VTwistSurface(name)
77 { 90 {
78 fHandedness = handedness; // +z = +ve, -z 91 fHandedness = handedness; // +z = +ve, -z = -ve
79 fAxis[0] = kRho; // in local coo 92 fAxis[0] = kRho; // in local coordinate system
80 fAxis[1] = kPhi; 93 fAxis[1] = kPhi;
81 G4int i = (handedness < 0 ? 0 : 1); 94 G4int i = (handedness < 0 ? 0 : 1);
82 fAxisMin[0] = EndInnerRadius[i]; // Inner- 95 fAxisMin[0] = EndInnerRadius[i]; // Inner-hype radius at z=0
83 fAxisMax[0] = EndOuterRadius[i]; // Outer- 96 fAxisMax[0] = EndOuterRadius[i]; // Outer-hype radius at z=0
84 fAxisMin[1] = -0.5*DPhi; 97 fAxisMin[1] = -0.5*DPhi;
85 fAxisMax[1] = -fAxisMin[1]; 98 fAxisMax[1] = -fAxisMin[1];
86 fCurrentNormal.normal.set(0, 0, (fHandednes 99 fCurrentNormal.normal.set(0, 0, (fHandedness < 0 ? -1 : 1));
87 // Unit vector, in local coordinate s 100 // Unit vector, in local coordinate system
88 fRot.rotateZ(EndPhi[i]); 101 fRot.rotateZ(EndPhi[i]);
89 fTrans.set(0, 0, EndZ[i]); 102 fTrans.set(0, 0, EndZ[i]);
90 fIsValidNorm = true; 103 fIsValidNorm = true;
91 104
92 SetCorners(); 105 SetCorners();
93 SetBoundaries(); 106 SetBoundaries();
94 107
95 fSurfaceArea = 0.5*DPhi * (EndOuterRadius[ 108 fSurfaceArea = 0.5*DPhi * (EndOuterRadius[i]*EndOuterRadius[i]
96 - EndInnerRadius[ << 109 - EndInnerRadius[i]*EndInnerRadius[i] ) ;
>> 110
97 } 111 }
98 112
>> 113
99 //============================================ 114 //=====================================================================
100 //* Fake default constructor ----------------- 115 //* Fake default constructor ------------------------------------------
101 116
102 G4TwistTubsFlatSide::G4TwistTubsFlatSide( __vo 117 G4TwistTubsFlatSide::G4TwistTubsFlatSide( __void__& a )
103 : G4VTwistSurface(a) << 118 : G4VTwistSurface(a), fSurfaceArea(0.)
104 { 119 {
105 } 120 }
106 121
>> 122
107 //============================================ 123 //=====================================================================
108 //* destructor ------------------------------- 124 //* destructor --------------------------------------------------------
109 125
110 G4TwistTubsFlatSide::~G4TwistTubsFlatSide() = << 126 G4TwistTubsFlatSide::~G4TwistTubsFlatSide()
>> 127 {
>> 128 }
111 129
112 //============================================ 130 //=====================================================================
113 //* GetNormal -------------------------------- 131 //* GetNormal ---------------------------------------------------------
114 132
115 G4ThreeVector G4TwistTubsFlatSide::GetNormal(c << 133 G4ThreeVector G4TwistTubsFlatSide::GetNormal(const G4ThreeVector & /* xx */ ,
116 << 134 G4bool isGlobal)
117 { 135 {
118 if (isGlobal) << 136 if (isGlobal) {
119 { <<
120 return ComputeGlobalDirection(fCurrentNo 137 return ComputeGlobalDirection(fCurrentNormal.normal);
121 } << 138 } else {
122 else <<
123 { <<
124 return fCurrentNormal.normal; 139 return fCurrentNormal.normal;
125 } 140 }
126 } 141 }
127 142
128 //============================================ 143 //=====================================================================
129 //* DistanceToSurface(p, v) ------------------ 144 //* DistanceToSurface(p, v) -------------------------------------------
130 145
131 G4int G4TwistTubsFlatSide::DistanceToSurface(c << 146 G4int G4TwistTubsFlatSide::DistanceToSurface(const G4ThreeVector &gp,
132 c << 147 const G4ThreeVector &gv,
133 << 148 G4ThreeVector gxx[],
134 << 149 G4double distance[],
135 << 150 G4int areacode[],
136 << 151 G4bool isvalid[],
137 << 152 EValidate validate)
138 { 153 {
139 fCurStatWithV.ResetfDone(validate, &gp, &gv 154 fCurStatWithV.ResetfDone(validate, &gp, &gv);
140 155
141 if (fCurStatWithV.IsDone()) << 156 if (fCurStatWithV.IsDone()) {
142 { << 157 G4int i;
143 for (G4int i=0; i<fCurStatWithV.GetNXX() << 158 for (i=0; i<fCurStatWithV.GetNXX(); i++) {
144 { <<
145 gxx[i] = fCurStatWithV.GetXX(i); 159 gxx[i] = fCurStatWithV.GetXX(i);
146 distance[i] = fCurStatWithV.GetDistan 160 distance[i] = fCurStatWithV.GetDistance(i);
147 areacode[i] = fCurStatWithV.GetAreaco 161 areacode[i] = fCurStatWithV.GetAreacode(i);
148 isvalid[i] = fCurStatWithV.IsValid(i 162 isvalid[i] = fCurStatWithV.IsValid(i);
149 } 163 }
150 return fCurStatWithV.GetNXX(); 164 return fCurStatWithV.GetNXX();
151 } << 165 } else {
152 else // initialize << 166 // initialize
153 { << 167 G4int i;
154 for (G4int i=0; i<2; ++i) << 168 for (i=0; i<2; i++) {
155 { <<
156 distance[i] = kInfinity; 169 distance[i] = kInfinity;
157 areacode[i] = sOutside; 170 areacode[i] = sOutside;
158 isvalid[i] = false; 171 isvalid[i] = false;
159 gxx[i].set(kInfinity, kInfinity, kInf 172 gxx[i].set(kInfinity, kInfinity, kInfinity);
160 } 173 }
161 } 174 }
162 175
163 G4ThreeVector p = ComputeLocalPoint(gp); 176 G4ThreeVector p = ComputeLocalPoint(gp);
164 G4ThreeVector v = ComputeLocalDirection(gv) 177 G4ThreeVector v = ComputeLocalDirection(gv);
165 178
166 // 179 //
167 // special case! 180 // special case!
168 // if p is on surface, distance = 0. 181 // if p is on surface, distance = 0.
169 // 182 //
170 183
171 if (std::fabs(p.z()) == 0.) // if p is o << 184 if (std::fabs(p.z()) == 0.) { // if p is on the plane
172 { <<
173 distance[0] = 0; 185 distance[0] = 0;
174 G4ThreeVector xx = p; 186 G4ThreeVector xx = p;
175 gxx[0] = ComputeGlobalPoint(xx); 187 gxx[0] = ComputeGlobalPoint(xx);
176 188
177 if (validate == kValidateWithTol) << 189 if (validate == kValidateWithTol) {
178 { <<
179 areacode[0] = GetAreaCode(xx); 190 areacode[0] = GetAreaCode(xx);
180 if (!IsOutside(areacode[0])) << 191 if (!IsOutside(areacode[0])) {
181 { <<
182 isvalid[0] = true; 192 isvalid[0] = true;
183 } 193 }
184 } << 194 } else if (validate == kValidateWithoutTol) {
185 else if (validate == kValidateWithoutTol <<
186 { <<
187 areacode[0] = GetAreaCode(xx, false); 195 areacode[0] = GetAreaCode(xx, false);
188 if (IsInside(areacode[0])) << 196 if (IsInside(areacode[0])) {
189 { <<
190 isvalid[0] = true; 197 isvalid[0] = true;
191 } 198 }
192 } << 199 } else { // kDontValidate
193 else // kDontValidate <<
194 { <<
195 areacode[0] = sInside; 200 areacode[0] = sInside;
196 isvalid[0] = true; 201 isvalid[0] = true;
197 } 202 }
>> 203
198 return 1; 204 return 1;
199 } 205 }
200 // 206 //
201 // special case end 207 // special case end
202 // 208 //
203 209
204 if (v.z() == 0) << 210 if (v.z() == 0) {
205 { << 211
206 fCurStatWithV.SetCurrentStatus(0, gxx[0] 212 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
207 isvalid[0 213 isvalid[0], 0, validate, &gp, &gv);
208 return 0; 214 return 0;
209 } 215 }
210 216
211 distance[0] = - (p.z() / v.z()); 217 distance[0] = - (p.z() / v.z());
212 218
213 G4ThreeVector xx = p + distance[0]*v; 219 G4ThreeVector xx = p + distance[0]*v;
214 gxx[0] = ComputeGlobalPoint(xx); 220 gxx[0] = ComputeGlobalPoint(xx);
215 221
216 if (validate == kValidateWithTol) << 222 if (validate == kValidateWithTol) {
217 { <<
218 areacode[0] = GetAreaCode(xx); 223 areacode[0] = GetAreaCode(xx);
219 if (!IsOutside(areacode[0])) << 224 if (!IsOutside(areacode[0])) {
220 { <<
221 if (distance[0] >= 0) isvalid[0] = tr 225 if (distance[0] >= 0) isvalid[0] = true;
222 } 226 }
223 } << 227 } else if (validate == kValidateWithoutTol) {
224 else if (validate == kValidateWithoutTol) <<
225 { <<
226 areacode[0] = GetAreaCode(xx, false); 228 areacode[0] = GetAreaCode(xx, false);
227 if (IsInside(areacode[0])) << 229 if (IsInside(areacode[0])) {
228 { <<
229 if (distance[0] >= 0) isvalid[0] = tr 230 if (distance[0] >= 0) isvalid[0] = true;
230 } 231 }
231 } << 232 } else { // kDontValidate
232 else // kDontValidate <<
233 { <<
234 areacode[0] = sInside; 233 areacode[0] = sInside;
235 if (distance[0] >= 0) isvalid[0] = tr 234 if (distance[0] >= 0) isvalid[0] = true;
236 } 235 }
237 236
238 fCurStatWithV.SetCurrentStatus(0, gxx[0], d 237 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
239 isvalid[0], 238 isvalid[0], 1, validate, &gp, &gv);
240 239
241 #ifdef G4TWISTDEBUG 240 #ifdef G4TWISTDEBUG
242 G4cerr << "ERROR - G4TwistTubsFlatSide::Dis 241 G4cerr << "ERROR - G4TwistTubsFlatSide::DistanceToSurface(p,v)" << G4endl;
243 G4cerr << " Name : " << GetNa 242 G4cerr << " Name : " << GetName() << G4endl;
244 G4cerr << " xx : " << xx << 243 G4cerr << " xx : " << xx << G4endl;
245 G4cerr << " gxx[0] : " << gxx[0 244 G4cerr << " gxx[0] : " << gxx[0] << G4endl;
246 G4cerr << " dist[0] : " << dista 245 G4cerr << " dist[0] : " << distance[0] << G4endl;
247 G4cerr << " areacode[0] : " << areac 246 G4cerr << " areacode[0] : " << areacode[0] << G4endl;
248 G4cerr << " isvalid[0] : " << isval 247 G4cerr << " isvalid[0] : " << isvalid[0] << G4endl;
249 #endif 248 #endif
250 return 1; 249 return 1;
251 } 250 }
252 251
253 //============================================ 252 //=====================================================================
254 //* DistanceToSurface(p) --------------------- 253 //* DistanceToSurface(p) ----------------------------------------------
255 254
256 G4int G4TwistTubsFlatSide::DistanceToSurface(c << 255 G4int G4TwistTubsFlatSide::DistanceToSurface(const G4ThreeVector &gp,
257 G 256 G4ThreeVector gxx[],
258 G 257 G4double distance[],
259 G 258 G4int areacode[])
260 { 259 {
261 // Calculate distance to plane in local coo 260 // Calculate distance to plane in local coordinate,
262 // then return distance and global intersec 261 // then return distance and global intersection points.
263 // 262 //
264 263
265 fCurStat.ResetfDone(kDontValidate, &gp); 264 fCurStat.ResetfDone(kDontValidate, &gp);
266 265
267 if (fCurStat.IsDone()) << 266 if (fCurStat.IsDone()) {
268 { << 267 G4int i;
269 for (G4int i=0; i<fCurStat.GetNXX(); ++i << 268 for (i=0; i<fCurStat.GetNXX(); i++) {
270 { <<
271 gxx[i] = fCurStat.GetXX(i); 269 gxx[i] = fCurStat.GetXX(i);
272 distance[i] = fCurStat.GetDistance(i) 270 distance[i] = fCurStat.GetDistance(i);
273 areacode[i] = fCurStat.GetAreacode(i) 271 areacode[i] = fCurStat.GetAreacode(i);
274 } 272 }
275 return fCurStat.GetNXX(); 273 return fCurStat.GetNXX();
276 } << 274 } else {
277 else // initialize << 275 // initialize
278 { << 276 G4int i;
279 for (auto i=0; i<2; ++i) << 277 for (i=0; i<2; i++) {
280 { <<
281 distance[i] = kInfinity; 278 distance[i] = kInfinity;
282 areacode[i] = sOutside; 279 areacode[i] = sOutside;
283 gxx[i].set(kInfinity, kInfinity, kInf 280 gxx[i].set(kInfinity, kInfinity, kInfinity);
284 } 281 }
285 } 282 }
286 283
287 G4ThreeVector p = ComputeLocalPoint(gp); 284 G4ThreeVector p = ComputeLocalPoint(gp);
288 G4ThreeVector xx; 285 G4ThreeVector xx;
289 286
290 // The plane is placed on origin with makin 287 // The plane is placed on origin with making its normal
291 // parallel to z-axis. 288 // parallel to z-axis.
292 if (std::fabs(p.z()) <= 0.5 * kCarTolerance << 289 if (std::fabs(p.z()) <= 0.5 * kCarTolerance) { // if p is on the plane, return 1
293 { << 290 distance[0] = 0;
294 distance[0] = 0; <<
295 xx = p; 291 xx = p;
296 } << 292 } else {
297 else <<
298 { <<
299 distance[0] = std::fabs(p.z()); 293 distance[0] = std::fabs(p.z());
300 xx.set(p.x(), p.y(), 0); 294 xx.set(p.x(), p.y(), 0);
301 } 295 }
302 296
303 gxx[0] = ComputeGlobalPoint(xx); 297 gxx[0] = ComputeGlobalPoint(xx);
304 areacode[0] = sInside; 298 areacode[0] = sInside;
305 G4bool isvalid = true; 299 G4bool isvalid = true;
306 fCurStat.SetCurrentStatus(0, gxx[0], distan 300 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
307 isvalid, 1, kDont 301 isvalid, 1, kDontValidate, &gp);
308 return 1; 302 return 1;
>> 303
309 } 304 }
310 305
311 //============================================ 306 //=====================================================================
312 //* GetAreaCode ------------------------------ 307 //* GetAreaCode -------------------------------------------------------
313 308
314 G4int G4TwistTubsFlatSide::GetAreaCode(const G 309 G4int G4TwistTubsFlatSide::GetAreaCode(const G4ThreeVector &xx,
315 G << 310 G4bool withTol)
316 { 311 {
317 const G4double rtol 312 const G4double rtol
318 = 0.5*G4GeometryTolerance::GetInstance()- 313 = 0.5*G4GeometryTolerance::GetInstance()->GetRadialTolerance();
319 314
320 G4int areacode = sInside; 315 G4int areacode = sInside;
321 316
322 if (fAxis[0] == kRho && fAxis[1] == kPhi) << 317 if (fAxis[0] == kRho && fAxis[1] == kPhi) {
323 { <<
324 G4int rhoaxis = 0; 318 G4int rhoaxis = 0;
>> 319 // G4int phiaxis = 0;
325 320
326 G4ThreeVector dphimin; // direction of 321 G4ThreeVector dphimin; // direction of phi-minimum boundary
327 G4ThreeVector dphimax; // direction of 322 G4ThreeVector dphimax; // direction of phi-maximum boundary
328 dphimin = GetCorner(sC0Max1Min); 323 dphimin = GetCorner(sC0Max1Min);
329 dphimax = GetCorner(sC0Max1Max); 324 dphimax = GetCorner(sC0Max1Max);
330 325
331 if (withTol) << 326 if (withTol) {
332 { << 327
333 G4bool isoutside = false; 328 G4bool isoutside = false;
334 329
335 // test boundary of rho-axis 330 // test boundary of rho-axis
336 331
337 if (xx.getRho() <= fAxisMin[rhoaxis] << 332 if (xx.getRho() <= fAxisMin[rhoaxis] + rtol) {
338 { << 333
339 areacode |= (sAxis0 & (sAxisRho|sA << 334 areacode |= (sAxis0 & (sAxisRho | sAxisMin)) | sBoundary; // rho-min
340 if (xx.getRho() < fAxisMin[rhoaxis 335 if (xx.getRho() < fAxisMin[rhoaxis] - rtol) isoutside = true;
341 336
342 } << 337 } else if (xx.getRho() >= fAxisMax[rhoaxis] - rtol) {
343 else if (xx.getRho() >= fAxisMax[rhoa << 338
344 { << 339 areacode |= (sAxis0 & (sAxisRho | sAxisMax)) | sBoundary; // rho-max
345 areacode |= (sAxis0 & (sAxisRho|sA << 340 if (xx.getRho() > fAxisMax[rhoaxis] + rtol) isoutside = true;
346 if (xx.getRho() > fAxisMax[rhoaxis << 341
347 } 342 }
348 343
349 // test boundary of phi-axis 344 // test boundary of phi-axis
350 345
351 if (AmIOnLeftSide(xx, dphimin) >= 0) << 346 if (AmIOnLeftSide(xx, dphimin) >= 0) { // xx is on dphimin
352 { << 347
353 areacode |= (sAxis1 & (sAxisPhi | 348 areacode |= (sAxis1 & (sAxisPhi | sAxisMin));
354 if ((areacode & sBoundary) != 0) << 349 if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner.
355 else areaco 350 else areacode |= sBoundary;
356 351
357 if (AmIOnLeftSide(xx, dphimin) > 0 352 if (AmIOnLeftSide(xx, dphimin) > 0) isoutside = true;
358 353
359 } << 354 } else if (AmIOnLeftSide(xx, dphimax) <= 0) { // xx is on dphimax
360 else if (AmIOnLeftSide(xx, dphimax) < << 355
361 { <<
362 areacode |= (sAxis1 & (sAxisPhi | 356 areacode |= (sAxis1 & (sAxisPhi | sAxisMax));
363 if ((areacode & sBoundary) != 0) << 357 if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner.
364 else areaco 358 else areacode |= sBoundary;
365 359
366 if (AmIOnLeftSide(xx, dphimax) < 0 << 360 if (AmIOnLeftSide(xx, dphimax) < 0) isoutside = true;
>> 361
367 } 362 }
368 363
369 // if isoutside = true, clear inside 364 // if isoutside = true, clear inside bit.
370 // if not on boundary, add axis infor 365 // if not on boundary, add axis information.
371 366
372 if (isoutside) << 367 if (isoutside) {
373 { <<
374 G4int tmpareacode = areacode & (~s 368 G4int tmpareacode = areacode & (~sInside);
375 areacode = tmpareacode; 369 areacode = tmpareacode;
376 } << 370 } else if ((areacode & sBoundary) != sBoundary) {
377 else if ((areacode & sBoundary) != sB <<
378 { <<
379 areacode |= (sAxis0 & sAxisRho) | 371 areacode |= (sAxis0 & sAxisRho) | (sAxis1 & sAxisPhi);
380 } 372 }
381 373
382 } << 374 } else {
383 else << 375
384 { <<
385 // out of boundary of rho-axis 376 // out of boundary of rho-axis
386 377
387 if (xx.getRho() < fAxisMin[rhoaxis]) << 378 if (xx.getRho() < fAxisMin[rhoaxis]) {
388 { <<
389 areacode |= (sAxis0 & (sAxisRho | 379 areacode |= (sAxis0 & (sAxisRho | sAxisMin)) | sBoundary;
390 } << 380 } else if (xx.getRho() > fAxisMax[rhoaxis]) {
391 else if (xx.getRho() > fAxisMax[rhoax <<
392 { <<
393 areacode |= (sAxis0 & (sAxisRho | 381 areacode |= (sAxis0 & (sAxisRho | sAxisMax)) | sBoundary;
394 } 382 }
395 383
396 // out of boundary of phi-axis 384 // out of boundary of phi-axis
397 385
398 if (AmIOnLeftSide(xx, dphimin, false) << 386 if (AmIOnLeftSide(xx, dphimin, false) >= 0) { // xx is leftside or
399 { << 387 areacode |= (sAxis1 & (sAxisPhi | sAxisMin)) ; // boundary of dphimin
400 areacode |= (sAxis1 & (sAxisPhi | s << 388 if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner.
401 if ((areacode & sBoundary) != 0) << 389 else areacode |= sBoundary;
402 else areacod <<
403 390
404 } << 391 } else if (AmIOnLeftSide(xx, dphimax, false) <= 0) { // xx is rightside or
405 else if (AmIOnLeftSide(xx, dphimax, f << 392 areacode |= (sAxis1 & (sAxisPhi | sAxisMax)) ; // boundary of dphimax
406 { << 393 if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner.
407 areacode |= (sAxis1 & (sAxisPhi | s << 394 else areacode |= sBoundary;
408 if ((areacode & sBoundary) != 0) <<
409 else areacod <<
410 395
411 } 396 }
412 397
413 if ((areacode & sBoundary) != sBounda 398 if ((areacode & sBoundary) != sBoundary) {
414 areacode |= (sAxis0 & sAxisRho) | 399 areacode |= (sAxis0 & sAxisRho) | (sAxis1 & sAxisPhi);
415 } 400 }
416 401
417 } 402 }
418 return areacode; 403 return areacode;
419 } << 404 } else {
420 else << 405
421 { <<
422 std::ostringstream message; 406 std::ostringstream message;
423 message << "Feature NOT implemented !" < 407 message << "Feature NOT implemented !" << G4endl
424 << " fAxis[0] = " << fAxi 408 << " fAxis[0] = " << fAxis[0] << G4endl
425 << " fAxis[1] = " << fAxi 409 << " fAxis[1] = " << fAxis[1];
426 G4Exception("G4TwistTubsFlatSide::GetAre 410 G4Exception("G4TwistTubsFlatSide::GetAreaCode()", "GeomSolids0001",
427 FatalException, message); 411 FatalException, message);
428 } 412 }
429 return areacode; 413 return areacode;
430 } 414 }
431 << 415
>> 416
432 //============================================ 417 //=====================================================================
433 //* SetCorners ------------------------------- 418 //* SetCorners --------------------------------------------------------
434 419
435 void G4TwistTubsFlatSide::SetCorners() 420 void G4TwistTubsFlatSide::SetCorners()
436 { 421 {
437 // Set Corner points in local coodinate. 422 // Set Corner points in local coodinate.
438 423
439 if (fAxis[0] == kRho && fAxis[1] == kPhi) << 424 if (fAxis[0] == kRho && fAxis[1] == kPhi) {
440 { << 425
441 G4int rhoaxis = 0; // kRho 426 G4int rhoaxis = 0; // kRho
442 G4int phiaxis = 1; // kPhi 427 G4int phiaxis = 1; // kPhi
443 428
444 G4double x, y, z; 429 G4double x, y, z;
445 // corner of Axis0min and Axis1min 430 // corner of Axis0min and Axis1min
446 x = fAxisMin[rhoaxis]*std::cos(fAxisM 431 x = fAxisMin[rhoaxis]*std::cos(fAxisMin[phiaxis]);
447 y = fAxisMin[rhoaxis]*std::sin(fAxisM 432 y = fAxisMin[rhoaxis]*std::sin(fAxisMin[phiaxis]);
448 z = 0; 433 z = 0;
449 SetCorner(sC0Min1Min, x, y, z); 434 SetCorner(sC0Min1Min, x, y, z);
450 // corner of Axis0max and Axis1min 435 // corner of Axis0max and Axis1min
451 x = fAxisMax[rhoaxis]*std::cos(fAxisM 436 x = fAxisMax[rhoaxis]*std::cos(fAxisMin[phiaxis]);
452 y = fAxisMax[rhoaxis]*std::sin(fAxisM 437 y = fAxisMax[rhoaxis]*std::sin(fAxisMin[phiaxis]);
453 z = 0; 438 z = 0;
454 SetCorner(sC0Max1Min, x, y, z); 439 SetCorner(sC0Max1Min, x, y, z);
455 // corner of Axis0max and Axis1max 440 // corner of Axis0max and Axis1max
456 x = fAxisMax[rhoaxis]*std::cos(fAxisM 441 x = fAxisMax[rhoaxis]*std::cos(fAxisMax[phiaxis]);
457 y = fAxisMax[rhoaxis]*std::sin(fAxisM 442 y = fAxisMax[rhoaxis]*std::sin(fAxisMax[phiaxis]);
458 z = 0; 443 z = 0;
459 SetCorner(sC0Max1Max, x, y, z); 444 SetCorner(sC0Max1Max, x, y, z);
460 // corner of Axis0min and Axis1max 445 // corner of Axis0min and Axis1max
461 x = fAxisMin[rhoaxis]*std::cos(fAxisM 446 x = fAxisMin[rhoaxis]*std::cos(fAxisMax[phiaxis]);
462 y = fAxisMin[rhoaxis]*std::sin(fAxisM 447 y = fAxisMin[rhoaxis]*std::sin(fAxisMax[phiaxis]);
463 z = 0; 448 z = 0;
464 SetCorner(sC0Min1Max, x, y, z); 449 SetCorner(sC0Min1Max, x, y, z);
465 450
466 } << 451 } else {
467 else <<
468 { <<
469 std::ostringstream message; 452 std::ostringstream message;
470 message << "Feature NOT implemented !" < 453 message << "Feature NOT implemented !" << G4endl
471 << " fAxis[0] = " << fAxi 454 << " fAxis[0] = " << fAxis[0] << G4endl
472 << " fAxis[1] = " << fAxi 455 << " fAxis[1] = " << fAxis[1];
473 G4Exception("G4TwistTubsFlatSide::SetCor 456 G4Exception("G4TwistTubsFlatSide::SetCorners()", "GeomSolids0001",
474 FatalException, message); 457 FatalException, message);
475 } 458 }
476 } 459 }
477 460
478 //============================================ 461 //=====================================================================
479 //* SetBoundaries() -------------------------- 462 //* SetBoundaries() ---------------------------------------------------
480 463
481 void G4TwistTubsFlatSide::SetBoundaries() 464 void G4TwistTubsFlatSide::SetBoundaries()
482 { 465 {
483 // Set direction-unit vector of phi-boundar 466 // Set direction-unit vector of phi-boundary-lines in local coodinate.
484 // Don't call the function twice. 467 // Don't call the function twice.
485 468
486 if (fAxis[0] == kRho && fAxis[1] == kPhi) << 469 if (fAxis[0] == kRho && fAxis[1] == kPhi) {
487 { << 470
488 G4ThreeVector direction; 471 G4ThreeVector direction;
489 // sAxis0 & sAxisMin 472 // sAxis0 & sAxisMin
490 direction = GetCorner(sC0Min1Max) - GetC 473 direction = GetCorner(sC0Min1Max) - GetCorner(sC0Min1Min);
491 direction = direction.unit(); 474 direction = direction.unit();
492 SetBoundary(sAxis0 & (sAxisPhi | sAxisMi 475 SetBoundary(sAxis0 & (sAxisPhi | sAxisMin), direction,
493 GetCorner(sC0Min1Min), sAxis 476 GetCorner(sC0Min1Min), sAxisPhi);
494 477
495 // sAxis0 & sAxisMax 478 // sAxis0 & sAxisMax
496 direction = GetCorner(sC0Max1Max) - GetC 479 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min);
497 direction = direction.unit(); 480 direction = direction.unit();
498 SetBoundary(sAxis0 & (sAxisPhi | sAxisMa 481 SetBoundary(sAxis0 & (sAxisPhi | sAxisMax), direction,
499 GetCorner(sC0Max1Min), sAxis 482 GetCorner(sC0Max1Min), sAxisPhi);
500 483
501 // sAxis1 & sAxisMin 484 // sAxis1 & sAxisMin
502 direction = GetCorner(sC0Max1Min) - GetC 485 direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
503 direction = direction.unit(); 486 direction = direction.unit();
504 SetBoundary(sAxis1 & (sAxisRho | sAxisMi 487 SetBoundary(sAxis1 & (sAxisRho | sAxisMin), direction,
505 GetCorner(sC0Min1Min), sAxis 488 GetCorner(sC0Min1Min), sAxisRho);
506 489
507 // sAxis1 & sAxisMax 490 // sAxis1 & sAxisMax
508 direction = GetCorner(sC0Max1Max) - GetC 491 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max);
509 direction = direction.unit(); 492 direction = direction.unit();
510 SetBoundary(sAxis1 & (sAxisRho | sAxisMa 493 SetBoundary(sAxis1 & (sAxisRho | sAxisMax), direction,
511 GetCorner(sC0Min1Max), sAxis 494 GetCorner(sC0Min1Max), sAxisPhi);
512 } << 495 } else {
513 else <<
514 { <<
515 std::ostringstream message; 496 std::ostringstream message;
516 message << "Feature NOT implemented !" < 497 message << "Feature NOT implemented !" << G4endl
517 << " fAxis[0] = " << fAxi 498 << " fAxis[0] = " << fAxis[0] << G4endl
518 << " fAxis[1] = " << fAxi 499 << " fAxis[1] = " << fAxis[1];
519 G4Exception("G4TwistTubsFlatSide::SetBou 500 G4Exception("G4TwistTubsFlatSide::SetBoundaries()", "GeomSolids0001",
520 FatalException, message); 501 FatalException, message);
521 } 502 }
522 } 503 }
523 504
524 //============================================ 505 //=====================================================================
525 //* GetFacets() ------------------------------ 506 //* GetFacets() -------------------------------------------------------
526 507
527 void G4TwistTubsFlatSide::GetFacets( G4int k, 508 void G4TwistTubsFlatSide::GetFacets( G4int k, G4int n, G4double xyz[][3],
528 G4int fac 509 G4int faces[][4], G4int iside )
529 { 510 {
>> 511
530 G4ThreeVector p ; 512 G4ThreeVector p ;
531 513
532 G4double rmin = fAxisMin[0] ; 514 G4double rmin = fAxisMin[0] ;
533 G4double rmax = fAxisMax[0] ; 515 G4double rmax = fAxisMax[0] ;
534 G4double phimin, phimax ; 516 G4double phimin, phimax ;
535 517
536 G4double r,phi ; 518 G4double r,phi ;
>> 519
>> 520 G4int i,j ;
>> 521
537 G4int nnode,nface ; 522 G4int nnode,nface ;
538 523
539 for ( G4int i = 0 ; i<n ; ++i ) << 524 for ( i = 0 ; i<n ; i++ ) {
540 { << 525
541 r = rmin + i*(rmax-rmin)/(n-1) ; 526 r = rmin + i*(rmax-rmin)/(n-1) ;
542 527
543 phimin = GetBoundaryMin(r) ; 528 phimin = GetBoundaryMin(r) ;
544 phimax = GetBoundaryMax(r) ; 529 phimax = GetBoundaryMax(r) ;
545 530
546 for ( G4int j = 0 ; j<k ; ++j ) << 531 for ( j = 0 ; j<k ; j++ )
547 { 532 {
548 phi = phimin + j*(phimax-phimin)/(k-1) ; 533 phi = phimin + j*(phimax-phimin)/(k-1) ;
549 534
550 nnode = GetNode(i,j,k,n,iside) ; 535 nnode = GetNode(i,j,k,n,iside) ;
551 p = SurfacePoint(phi,r,true) ; // surfa 536 p = SurfacePoint(phi,r,true) ; // surface point in global coord.system
552 537
553 xyz[nnode][0] = p.x() ; 538 xyz[nnode][0] = p.x() ;
554 xyz[nnode][1] = p.y() ; 539 xyz[nnode][1] = p.y() ;
555 xyz[nnode][2] = p.z() ; 540 xyz[nnode][2] = p.z() ;
556 541
557 if ( i<n-1 && j<k-1 ) // conterclock << 542 if ( i<n-1 && j<k-1 ) { // conterclock wise filling
558 { << 543
559 nface = GetFace(i,j,k,n,iside) ; 544 nface = GetFace(i,j,k,n,iside) ;
560 545
561 if (fHandedness < 0) // lower side << 546 if (fHandedness < 0) { // lower side
562 { << 547 faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,-1) * ( GetNode(i ,j ,k,n,iside)+1) ;
563 faces[nface][0] = GetEdgeVisibility( << 548 faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,-1) * ( GetNode(i ,j+1,k,n,iside)+1) ;
564 * ( GetNode(i ,j , << 549 faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,-1) * ( GetNode(i+1,j+1,k,n,iside)+1) ;
565 faces[nface][1] = GetEdgeVisibility( << 550 faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,-1) * ( GetNode(i+1,j ,k,n,iside)+1) ;
566 * ( GetNode(i ,j+1, << 551 } else { // upper side
567 faces[nface][2] = GetEdgeVisibility( << 552 faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,1) * ( GetNode(i ,j ,k,n,iside)+1) ;
568 * ( GetNode(i+1,j+1, << 553 faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,1) * ( GetNode(i+1,j ,k,n,iside)+1) ;
569 faces[nface][3] = GetEdgeVisibility( << 554 faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,1) * ( GetNode(i+1,j+1,k,n,iside)+1) ;
570 * ( GetNode(i+1,j , << 555 faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,1) * ( GetNode(i ,j+1,k,n,iside)+1) ;
571 } << 556
572 else // upper side <<
573 { <<
574 faces[nface][0] = GetEdgeVisibility( <<
575 * ( GetNode(i ,j , <<
576 faces[nface][1] = GetEdgeVisibility( <<
577 * ( GetNode(i+1,j , <<
578 faces[nface][2] = GetEdgeVisibility( <<
579 * ( GetNode(i+1,j+1, <<
580 faces[nface][3] = GetEdgeVisibility( <<
581 * ( GetNode(i ,j+1, <<
582 } 557 }
>> 558
>> 559
>> 560
583 } 561 }
584 } 562 }
585 } 563 }
586 } 564 }
587 565