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