| 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,v 1.7 2007/05/23 09:31:02 gcosmo Exp $
29 // 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), << 28 // GEANT4 tag $Name: geant4-09-01-patch-03 $
30 // from original version in Jupi << 29 //
>> 30 //
>> 31 // --------------------------------------------------------------------
>> 32 // GEANT 4 class source file
>> 33 //
>> 34 //
>> 35 // G4TwistTubsFlatSide.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 44
33 #include "G4TwistTubsFlatSide.hh" 45 #include "G4TwistTubsFlatSide.hh"
34 #include "G4GeometryTolerance.hh" 46 #include "G4GeometryTolerance.hh"
35 47
36 //============================================ 48 //=====================================================================
37 //* constructors ----------------------------- 49 //* constructors ------------------------------------------------------
38 50
39 G4TwistTubsFlatSide::G4TwistTubsFlatSide(const << 51 G4TwistTubsFlatSide::G4TwistTubsFlatSide(const G4String &name,
40 const G4RotationM << 52 const G4RotationMatrix &rot,
41 const G4ThreeVect << 53 const G4ThreeVector &tlate,
42 const G4ThreeVect << 54 const G4ThreeVector &n,
43 const EAxis 55 const EAxis axis0 ,
44 const EAxis 56 const EAxis axis1 ,
45 G4double 57 G4double axis0min,
46 G4double 58 G4double axis1min,
47 G4double 59 G4double axis0max,
48 G4double 60 G4double axis1max )
49 : G4VTwistSurface(name, rot, tlate, 0, axis0 61 : G4VTwistSurface(name, rot, tlate, 0, axis0, axis1,
50 axis0min, axis1min, axis0m << 62 axis0min, axis1min, axis0max, axis1max)
51 { 63 {
52 if (axis0 == kPhi && axis1 == kRho) << 64 if (axis0 == kPhi && axis1 == kRho) {
53 { << 65 G4Exception("G4TwistTubsFlatSide::G4TwistTubsFlatSide()", "InvalidSetup",
54 G4Exception("G4TwistTubsFlatSide::G4Twis << 66 FatalException, "Should swap axis0 and axis1!");
55 "GeomSolids0002", FatalError <<
56 "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)
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
>> 313 static const G4double rtol
318 = 0.5*G4GeometryTolerance::GetInstance()- 314 = 0.5*G4GeometryTolerance::GetInstance()->GetRadialTolerance();
319 315
320 G4int areacode = sInside; 316 G4int areacode = sInside;
321 317
322 if (fAxis[0] == kRho && fAxis[1] == kPhi) << 318 if (fAxis[0] == kRho && fAxis[1] == kPhi) {
323 { <<
324 G4int rhoaxis = 0; 319 G4int rhoaxis = 0;
>> 320 // G4int phiaxis = 0;
325 321
326 G4ThreeVector dphimin; // direction of 322 G4ThreeVector dphimin; // direction of phi-minimum boundary
327 G4ThreeVector dphimax; // direction of 323 G4ThreeVector dphimax; // direction of phi-maximum boundary
328 dphimin = GetCorner(sC0Max1Min); 324 dphimin = GetCorner(sC0Max1Min);
329 dphimax = GetCorner(sC0Max1Max); 325 dphimax = GetCorner(sC0Max1Max);
330 326
331 if (withTol) << 327 if (withTol) {
332 { << 328
333 G4bool isoutside = false; 329 G4bool isoutside = false;
334 330
335 // test boundary of rho-axis 331 // test boundary of rho-axis
336 332
337 if (xx.getRho() <= fAxisMin[rhoaxis] << 333 if (xx.getRho() <= fAxisMin[rhoaxis] + rtol) {
338 { << 334
339 areacode |= (sAxis0 & (sAxisRho|sA << 335 areacode |= (sAxis0 & (sAxisRho | sAxisMin)) | sBoundary; // rho-min
340 if (xx.getRho() < fAxisMin[rhoaxis 336 if (xx.getRho() < fAxisMin[rhoaxis] - rtol) isoutside = true;
341 337
342 } << 338 } else if (xx.getRho() >= fAxisMax[rhoaxis] - rtol) {
343 else if (xx.getRho() >= fAxisMax[rhoa << 339
344 { << 340 areacode |= (sAxis0 & (sAxisRho | sAxisMax)) | sBoundary; // rho-max
345 areacode |= (sAxis0 & (sAxisRho|sA << 341 if (xx.getRho() > fAxisMax[rhoaxis] + rtol) isoutside = true;
346 if (xx.getRho() > fAxisMax[rhoaxis << 342
347 } 343 }
348 344
349 // test boundary of phi-axis 345 // test boundary of phi-axis
350 346
351 if (AmIOnLeftSide(xx, dphimin) >= 0) << 347 if (AmIOnLeftSide(xx, dphimin) >= 0) { // xx is on dphimin
352 { << 348
353 areacode |= (sAxis1 & (sAxisPhi | 349 areacode |= (sAxis1 & (sAxisPhi | sAxisMin));
354 if ((areacode & sBoundary) != 0) << 350 if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner.
355 else areaco 351 else areacode |= sBoundary;
356 352
357 if (AmIOnLeftSide(xx, dphimin) > 0 353 if (AmIOnLeftSide(xx, dphimin) > 0) isoutside = true;
358 354
359 } << 355 } else if (AmIOnLeftSide(xx, dphimax) <= 0) { // xx is on dphimax
360 else if (AmIOnLeftSide(xx, dphimax) < << 356
361 { <<
362 areacode |= (sAxis1 & (sAxisPhi | 357 areacode |= (sAxis1 & (sAxisPhi | sAxisMax));
363 if ((areacode & sBoundary) != 0) << 358 if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner.
364 else areaco 359 else areacode |= sBoundary;
365 360
366 if (AmIOnLeftSide(xx, dphimax) < 0 << 361 if (AmIOnLeftSide(xx, dphimax) < 0) isoutside = true;
>> 362
367 } 363 }
368 364
369 // if isoutside = true, clear inside 365 // if isoutside = true, clear inside bit.
370 // if not on boundary, add axis infor 366 // if not on boundary, add axis information.
371 367
372 if (isoutside) << 368 if (isoutside) {
373 { <<
374 G4int tmpareacode = areacode & (~s 369 G4int tmpareacode = areacode & (~sInside);
375 areacode = tmpareacode; 370 areacode = tmpareacode;
376 } << 371 } else if ((areacode & sBoundary) != sBoundary) {
377 else if ((areacode & sBoundary) != sB <<
378 { <<
379 areacode |= (sAxis0 & sAxisRho) | 372 areacode |= (sAxis0 & sAxisRho) | (sAxis1 & sAxisPhi);
380 } 373 }
381 374
382 } << 375 } else {
383 else << 376
384 { <<
385 // out of boundary of rho-axis 377 // out of boundary of rho-axis
386 378
387 if (xx.getRho() < fAxisMin[rhoaxis]) << 379 if (xx.getRho() < fAxisMin[rhoaxis]) {
388 { <<
389 areacode |= (sAxis0 & (sAxisRho | 380 areacode |= (sAxis0 & (sAxisRho | sAxisMin)) | sBoundary;
390 } << 381 } else if (xx.getRho() > fAxisMax[rhoaxis]) {
391 else if (xx.getRho() > fAxisMax[rhoax <<
392 { <<
393 areacode |= (sAxis0 & (sAxisRho | 382 areacode |= (sAxis0 & (sAxisRho | sAxisMax)) | sBoundary;
394 } 383 }
395 384
396 // out of boundary of phi-axis 385 // out of boundary of phi-axis
397 386
398 if (AmIOnLeftSide(xx, dphimin, false) << 387 if (AmIOnLeftSide(xx, dphimin, false) >= 0) { // xx is leftside or
399 { << 388 areacode |= (sAxis1 & (sAxisPhi | sAxisMin)) ; // boundary of dphimin
400 areacode |= (sAxis1 & (sAxisPhi | s << 389 if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner.
401 if ((areacode & sBoundary) != 0) << 390 else areacode |= sBoundary;
402 else areacod <<
403 391
404 } << 392 } else if (AmIOnLeftSide(xx, dphimax, false) <= 0) { // xx is rightside or
405 else if (AmIOnLeftSide(xx, dphimax, f << 393 areacode |= (sAxis1 & (sAxisPhi | sAxisMax)) ; // boundary of dphimax
406 { << 394 if (areacode & sBoundary) areacode |= sCorner; // xx is on the corner.
407 areacode |= (sAxis1 & (sAxisPhi | s << 395 else areacode |= sBoundary;
408 if ((areacode & sBoundary) != 0) <<
409 else areacod <<
410 396
411 } 397 }
412 398
413 if ((areacode & sBoundary) != sBounda 399 if ((areacode & sBoundary) != sBoundary) {
414 areacode |= (sAxis0 & sAxisRho) | 400 areacode |= (sAxis0 & sAxisRho) | (sAxis1 & sAxisPhi);
415 } 401 }
416 402
417 } 403 }
418 return areacode; 404 return areacode;
419 } << 405 } else {
420 else << 406
421 { << 407 G4cerr << "ERROR - G4TwistTubsFlatSide::GetAreaCode()" << G4endl
422 std::ostringstream message; << 408 << " fAxis[0] = " << fAxis[0] << G4endl
423 message << "Feature NOT implemented !" < << 409 << " fAxis[1] = " << fAxis[1] << G4endl;
424 << " fAxis[0] = " << fAxi << 410 G4Exception("G4TwistTubsFlatSide::GetAreaCode()", "NotImplemented",
425 << " fAxis[1] = " << fAxi << 411 FatalException, "Feature NOT implemented !");
426 G4Exception("G4TwistTubsFlatSide::GetAre <<
427 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 << 452 G4cerr << "ERROR - G4TwistTubsFlatSide::SetCorners()" << G4endl
468 { << 453 << " fAxis[0] = " << fAxis[0] << G4endl
469 std::ostringstream message; << 454 << " fAxis[1] = " << fAxis[1] << G4endl;
470 message << "Feature NOT implemented !" < << 455 G4Exception("G4TwistTubsFlatSide::SetCorners()", "NotImplemented",
471 << " fAxis[0] = " << fAxi << 456 FatalException, "Feature NOT implemented !");
472 << " fAxis[1] = " << fAxi <<
473 G4Exception("G4TwistTubsFlatSide::SetCor <<
474 FatalException, message); <<
475 } 457 }
476 } 458 }
477 459
478 //============================================ 460 //=====================================================================
479 //* SetBoundaries() -------------------------- 461 //* SetBoundaries() ---------------------------------------------------
480 462
481 void G4TwistTubsFlatSide::SetBoundaries() 463 void G4TwistTubsFlatSide::SetBoundaries()
482 { 464 {
483 // Set direction-unit vector of phi-boundar 465 // Set direction-unit vector of phi-boundary-lines in local coodinate.
484 // Don't call the function twice. 466 // Don't call the function twice.
485 467
486 if (fAxis[0] == kRho && fAxis[1] == kPhi) << 468 if (fAxis[0] == kRho && fAxis[1] == kPhi) {
487 { << 469
488 G4ThreeVector direction; 470 G4ThreeVector direction;
489 // sAxis0 & sAxisMin 471 // sAxis0 & sAxisMin
490 direction = GetCorner(sC0Min1Max) - GetC 472 direction = GetCorner(sC0Min1Max) - GetCorner(sC0Min1Min);
491 direction = direction.unit(); 473 direction = direction.unit();
492 SetBoundary(sAxis0 & (sAxisPhi | sAxisMi 474 SetBoundary(sAxis0 & (sAxisPhi | sAxisMin), direction,
493 GetCorner(sC0Min1Min), sAxis 475 GetCorner(sC0Min1Min), sAxisPhi);
494 476
495 // sAxis0 & sAxisMax 477 // sAxis0 & sAxisMax
496 direction = GetCorner(sC0Max1Max) - GetC 478 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min);
497 direction = direction.unit(); 479 direction = direction.unit();
498 SetBoundary(sAxis0 & (sAxisPhi | sAxisMa 480 SetBoundary(sAxis0 & (sAxisPhi | sAxisMax), direction,
499 GetCorner(sC0Max1Min), sAxis 481 GetCorner(sC0Max1Min), sAxisPhi);
500 482
501 // sAxis1 & sAxisMin 483 // sAxis1 & sAxisMin
502 direction = GetCorner(sC0Max1Min) - GetC 484 direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
503 direction = direction.unit(); 485 direction = direction.unit();
504 SetBoundary(sAxis1 & (sAxisRho | sAxisMi 486 SetBoundary(sAxis1 & (sAxisRho | sAxisMin), direction,
505 GetCorner(sC0Min1Min), sAxis 487 GetCorner(sC0Min1Min), sAxisRho);
506 488
507 // sAxis1 & sAxisMax 489 // sAxis1 & sAxisMax
508 direction = GetCorner(sC0Max1Max) - GetC 490 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max);
509 direction = direction.unit(); 491 direction = direction.unit();
510 SetBoundary(sAxis1 & (sAxisRho | sAxisMa 492 SetBoundary(sAxis1 & (sAxisRho | sAxisMax), direction,
511 GetCorner(sC0Min1Max), sAxis 493 GetCorner(sC0Min1Max), sAxisPhi);
512 } << 494 } else {
513 else << 495 G4cerr << "ERROR - G4TwistTubsFlatSide::SetBoundaries()" << G4endl
514 { << 496 << " fAxis[0] = " << fAxis[0] << G4endl
515 std::ostringstream message; << 497 << " fAxis[1] = " << fAxis[1] << G4endl;
516 message << "Feature NOT implemented !" < << 498 G4Exception("G4TwistTubsFlatSide::SetBoundaries()", "NotImplemented",
517 << " fAxis[0] = " << fAxi << 499 FatalException, "Feature NOT implemented !");
518 << " fAxis[1] = " << fAxi <<
519 G4Exception("G4TwistTubsFlatSide::SetBou <<
520 FatalException, message); <<
521 } 500 }
522 } 501 }
523 502
524 //============================================ 503 //=====================================================================
525 //* GetFacets() ------------------------------ 504 //* GetFacets() -------------------------------------------------------
526 505
527 void G4TwistTubsFlatSide::GetFacets( G4int k, << 506 void G4TwistTubsFlatSide::GetFacets( G4int m, G4int n, G4double xyz[][3],
528 G4int fac 507 G4int faces[][4], G4int iside )
529 { 508 {
>> 509
530 G4ThreeVector p ; 510 G4ThreeVector p ;
531 511
532 G4double rmin = fAxisMin[0] ; 512 G4double rmin = fAxisMin[0] ;
533 G4double rmax = fAxisMax[0] ; 513 G4double rmax = fAxisMax[0] ;
534 G4double phimin, phimax ; 514 G4double phimin, phimax ;
535 515
536 G4double r,phi ; 516 G4double r,phi ;
>> 517
>> 518 G4int i,j ;
>> 519
537 G4int nnode,nface ; 520 G4int nnode,nface ;
538 521
539 for ( G4int i = 0 ; i<n ; ++i ) << 522 for ( i = 0 ; i<n ; i++ ) {
540 { << 523
541 r = rmin + i*(rmax-rmin)/(n-1) ; 524 r = rmin + i*(rmax-rmin)/(n-1) ;
542 525
543 phimin = GetBoundaryMin(r) ; 526 phimin = GetBoundaryMin(r) ;
544 phimax = GetBoundaryMax(r) ; 527 phimax = GetBoundaryMax(r) ;
545 528
546 for ( G4int j = 0 ; j<k ; ++j ) << 529 for ( j = 0 ; j<m ; j++ )
547 { 530 {
548 phi = phimin + j*(phimax-phimin)/(k-1) ; << 531 phi = phimin + j*(phimax-phimin)/(m-1) ;
549 532
550 nnode = GetNode(i,j,k,n,iside) ; << 533 nnode = GetNode(i,j,m,n,iside) ;
551 p = SurfacePoint(phi,r,true) ; // surfa 534 p = SurfacePoint(phi,r,true) ; // surface point in global coord.system
552 535
553 xyz[nnode][0] = p.x() ; 536 xyz[nnode][0] = p.x() ;
554 xyz[nnode][1] = p.y() ; 537 xyz[nnode][1] = p.y() ;
555 xyz[nnode][2] = p.z() ; 538 xyz[nnode][2] = p.z() ;
556 539
557 if ( i<n-1 && j<k-1 ) // conterclock << 540 if ( i<n-1 && j<m-1 ) { // conterclock wise filling
558 { << 541
559 nface = GetFace(i,j,k,n,iside) ; << 542 nface = GetFace(i,j,m,n,iside) ;
560 << 543
561 if (fHandedness < 0) // lower side << 544 if (fHandedness < 0) { // lower side
562 { << 545 faces[nface][0] = GetEdgeVisibility(i,j,m,n,0,-1) * ( GetNode(i ,j ,m,n,iside)+1) ;
563 faces[nface][0] = GetEdgeVisibility( << 546 faces[nface][1] = GetEdgeVisibility(i,j,m,n,1,-1) * ( GetNode(i ,j+1,m,n,iside)+1) ;
564 * ( GetNode(i ,j , << 547 faces[nface][2] = GetEdgeVisibility(i,j,m,n,2,-1) * ( GetNode(i+1,j+1,m,n,iside)+1) ;
565 faces[nface][1] = GetEdgeVisibility( << 548 faces[nface][3] = GetEdgeVisibility(i,j,m,n,3,-1) * ( GetNode(i+1,j ,m,n,iside)+1) ;
566 * ( GetNode(i ,j+1, << 549 } else { // upper side
567 faces[nface][2] = GetEdgeVisibility( << 550 faces[nface][0] = GetEdgeVisibility(i,j,m,n,0,1) * ( GetNode(i ,j ,m,n,iside)+1) ;
568 * ( GetNode(i+1,j+1, << 551 faces[nface][1] = GetEdgeVisibility(i,j,m,n,1,1) * ( GetNode(i+1,j ,m,n,iside)+1) ;
569 faces[nface][3] = GetEdgeVisibility( << 552 faces[nface][2] = GetEdgeVisibility(i,j,m,n,2,1) * ( GetNode(i+1,j+1,m,n,iside)+1) ;
570 * ( GetNode(i+1,j , << 553 faces[nface][3] = GetEdgeVisibility(i,j,m,n,3,1) * ( GetNode(i ,j+1,m,n,iside)+1) ;
571 } << 554
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 } 555 }
>> 556
>> 557
>> 558
583 } 559 }
584 } 560 }
585 } 561 }
586 } 562 }
587 563