Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // G4TwistTrapFlatSide implementation 27 // 28 // Author: 30-Aug-2002 - O.Link (Oliver.Link@c 29 // ------------------------------------------- 30 31 #include "G4TwistTrapFlatSide.hh" 32 33 //============================================ 34 //* constructors ----------------------------- 35 36 G4TwistTrapFlatSide::G4TwistTrapFlatSide( cons 37 38 39 40 41 42 43 44 45 46 : G4VTwistSurface(name) 47 { 48 fHandedness = handedness; // +z = +ve, -z 49 50 fDx1 = pDx1 ; 51 fDx2 = pDx2 ; 52 fDy = pDy ; 53 fDz = pDz ; 54 fAlpha = pAlpha ; 55 fTAlph = std::tan(fAlpha) ; 56 fPhi = pPhi ; 57 fTheta = pTheta ; 58 59 fdeltaX = 2 * fDz * std::tan(fTheta) * std: 60 // dx in surface equation 61 fdeltaY = 2 * fDz * std::tan(fTheta) * std: 62 // dy in surface equation 63 64 fPhiTwist = PhiTwist ; 65 66 fCurrentNormal.normal.set( 0, 0, (fHandedne 67 // Unit vector, in local coordinate s 68 fRot.rotateZ( fHandedness > 0 69 ? 0.5 * fPhiTwist 70 : -0.5 * fPhiTwist ); 71 72 fTrans.set( 73 fHandedness > 0 ? 0.5*fdeltaX : 74 fHandedness > 0 ? 0.5*fdeltaY : 75 fHandedness > 0 ? fDz : -fDz ) ; 76 77 fIsValidNorm = true; 78 79 80 fAxis[0] = kXAxis ; 81 fAxis[1] = kYAxis ; 82 fAxisMin[0] = kInfinity ; // x-Axis cannot 83 fAxisMax[0] = kInfinity ; // depends on y 84 fAxisMin[1] = -fDy ; // y - axis 85 fAxisMax[1] = fDy ; 86 87 SetCorners(); 88 SetBoundaries(); 89 } 90 91 92 //============================================ 93 //* Fake default constructor ----------------- 94 95 G4TwistTrapFlatSide::G4TwistTrapFlatSide( __vo 96 : G4VTwistSurface(a) 97 { 98 } 99 100 101 //============================================ 102 //* destructor ------------------------------- 103 104 G4TwistTrapFlatSide::~G4TwistTrapFlatSide() = 105 106 //============================================ 107 //* GetNormal -------------------------------- 108 109 G4ThreeVector G4TwistTrapFlatSide::GetNormal(c 110 G 111 { 112 if (isGlobal) 113 { 114 return ComputeGlobalDirection(fCurrentNo 115 } 116 else 117 { 118 return fCurrentNormal.normal; 119 } 120 } 121 122 //============================================ 123 //* DistanceToSurface(p, v) ------------------ 124 125 G4int G4TwistTrapFlatSide::DistanceToSurface(c 126 c 127 128 129 130 131 132 { 133 fCurStatWithV.ResetfDone(validate, &gp, &gv 134 135 if (fCurStatWithV.IsDone()) 136 { 137 for (G4int i=0; i<fCurStatWithV.GetNXX() 138 { 139 gxx[i] = fCurStatWithV.GetXX(i); 140 distance[i] = fCurStatWithV.GetDistan 141 areacode[i] = fCurStatWithV.GetAreaco 142 isvalid[i] = fCurStatWithV.IsValid(i 143 } 144 return fCurStatWithV.GetNXX(); 145 } 146 else // initialize 147 { 148 for (auto i=0; i<2; ++i) 149 { 150 distance[i] = kInfinity; 151 areacode[i] = sOutside; 152 isvalid[i] = false; 153 gxx[i].set(kInfinity, kInfinity, kInf 154 } 155 } 156 157 G4ThreeVector p = ComputeLocalPoint(gp); 158 G4ThreeVector v = ComputeLocalDirection(gv) 159 160 // 161 // special case! 162 // if p is on surface, distance = 0. 163 // 164 165 if (std::fabs(p.z()) == 0.) // if p is 166 { 167 distance[0] = 0; 168 G4ThreeVector xx = p; 169 gxx[0] = ComputeGlobalPoint(xx); 170 171 if (validate == kValidateWithTol) 172 { 173 areacode[0] = GetAreaCode(xx); 174 if (!IsOutside(areacode[0])) 175 { 176 isvalid[0] = true; 177 } 178 } 179 else if (validate == kValidateWithoutTol 180 { 181 areacode[0] = GetAreaCode(xx, false); 182 if (IsInside(areacode[0])) 183 { 184 isvalid[0] = true; 185 } 186 } 187 else // kDontValidate 188 { 189 areacode[0] = sInside; 190 isvalid[0] = true; 191 } 192 return 1; 193 } 194 // 195 // special case end 196 // 197 198 if (v.z() == 0) { 199 200 fCurStatWithV.SetCurrentStatus(0, gxx[0] 201 isvalid[0 202 return 0; 203 } 204 205 distance[0] = - (p.z() / v.z()); 206 207 G4ThreeVector xx = p + distance[0]*v; 208 gxx[0] = ComputeGlobalPoint(xx); 209 210 if (validate == kValidateWithTol) 211 { 212 areacode[0] = GetAreaCode(xx); 213 if (!IsOutside(areacode[0])) 214 { 215 if (distance[0] >= 0) isvalid[0] = tr 216 } 217 } 218 else if (validate == kValidateWithoutTol) 219 { 220 areacode[0] = GetAreaCode(xx, false); 221 if (IsInside(areacode[0])) 222 { 223 if (distance[0] >= 0) isvalid[0] = tr 224 } 225 } 226 else // kDontValidate 227 { 228 areacode[0] = sInside; 229 if (distance[0] >= 0) isvalid[0] = tr 230 } 231 232 fCurStatWithV.SetCurrentStatus(0, gxx[0], d 233 isvalid[0], 234 235 #ifdef G4TWISTDEBUG 236 G4cerr << "ERROR - G4TwistTrapFlatSide::Dis 237 G4cerr << " Name : " << GetNa 238 G4cerr << " xx : " << xx << 239 G4cerr << " gxx[0] : " << gxx[0 240 G4cerr << " dist[0] : " << dista 241 G4cerr << " areacode[0] : " << areac 242 G4cerr << " isvalid[0] : " << isval 243 #endif 244 return 1; 245 } 246 247 //============================================ 248 //* DistanceToSurface(p) --------------------- 249 250 G4int G4TwistTrapFlatSide::DistanceToSurface(c 251 G 252 G 253 G 254 { 255 // Calculate distance to plane in local coo 256 // then return distance and global intersec 257 // 258 259 fCurStat.ResetfDone(kDontValidate, &gp); 260 261 if (fCurStat.IsDone()) 262 { 263 for (G4int i=0; i<fCurStat.GetNXX(); ++i 264 { 265 gxx[i] = fCurStat.GetXX(i); 266 distance[i] = fCurStat.GetDistance(i) 267 areacode[i] = fCurStat.GetAreacode(i) 268 } 269 return fCurStat.GetNXX(); 270 } 271 else // initialize 272 { 273 for (auto i=0; i<2; ++i) 274 { 275 distance[i] = kInfinity; 276 areacode[i] = sOutside; 277 gxx[i].set(kInfinity, kInfinity, kInf 278 } 279 } 280 281 G4ThreeVector p = ComputeLocalPoint(gp); 282 G4ThreeVector xx; 283 284 // The plane is placed on origin with makin 285 // parallel to z-axis. 286 if (std::fabs(p.z()) <= 0.5 * kCarTolerance 287 { // if p is on the plane, return 1 288 distance[0] = 0; 289 xx = p; 290 } 291 else 292 { 293 distance[0] = std::fabs(p.z()); 294 xx.set(p.x(), p.y(), 0); 295 } 296 297 gxx[0] = ComputeGlobalPoint(xx); 298 areacode[0] = sInside; 299 G4bool isvalid = true; 300 fCurStat.SetCurrentStatus(0, gxx[0], distan 301 isvalid, 1, kDont 302 return 1; 303 304 } 305 306 //============================================ 307 //* GetAreaCode() ---------------------------- 308 309 G4int G4TwistTrapFlatSide::GetAreaCode(const G 310 G 311 { 312 313 static const G4double ctol = 0.5 * kCarToler 314 G4int areacode = sInside; 315 316 if (fAxis[0] == kXAxis && fAxis[1] == kYAxis 317 { 318 G4int yaxis = 1; 319 320 G4double wmax = xAxisMax(xx.y(), fTAlph ) 321 G4double wmin = -xAxisMax(xx.y(), -fTAlph 322 323 if (withTol) 324 { 325 G4bool isoutside = false; 326 327 // test boundary of x-axis 328 329 if (xx.x() < wmin + ctol) 330 { 331 areacode |= (sAxis0 & (sAxisX | sAxisM 332 if (xx.x() <= wmin - ctol) isoutside = 333 334 } 335 else if (xx.x() > wmax - ctol) 336 { 337 areacode |= (sAxis0 & (sAxisX | sAxisM 338 if (xx.x() >= wmax + ctol) isoutside 339 } 340 341 // test boundary of y-axis 342 343 if (xx.y() < fAxisMin[yaxis] + ctol) 344 { 345 areacode |= (sAxis1 & (sAxisY | sAxisM 346 347 if ((areacode & sBoundary) != 0) are 348 else areacode | 349 if (xx.y() <= fAxisMin[yaxis] - ctol) 350 351 } 352 else if (xx.y() > fAxisMax[yaxis] - ctol 353 { 354 areacode |= (sAxis1 & (sAxisY | sAxisM 355 356 if ((areacode & sBoundary) != 0) are 357 else areacode | 358 if (xx.y() >= fAxisMax[yaxis] + ctol) 359 } 360 361 // if isoutside = true, clear inside bit 362 // if not on boundary, add axis informat 363 364 if (isoutside) 365 { 366 G4int tmpareacode = areacode & (~sInsi 367 areacode = tmpareacode; 368 } 369 else if ((areacode & sBoundary) != sBoun 370 { 371 areacode |= (sAxis0 & sAxisX) | (sAxis 372 } 373 } 374 else 375 { 376 // boundary of x-axis 377 378 if (xx.x() < wmin ) 379 { 380 areacode |= (sAxis0 & (sAxisX | sAxisM 381 } 382 else if (xx.x() > wmax) 383 { 384 areacode |= (sAxis0 & (sAxisX | sAxisM 385 } 386 387 // boundary of y-axis 388 389 if (xx.y() < fAxisMin[yaxis]) 390 { 391 areacode |= (sAxis1 & (sAxisY | sAxisM 392 if ((areacode & sBoundary) != 0) are 393 else areacode | 394 395 } 396 else if (xx.y() > fAxisMax[yaxis]) 397 { 398 areacode |= (sAxis1 & (sAxisY | sAxisM 399 if ((areacode & sBoundary) != 0) are 400 else areacode | 401 } 402 403 if ((areacode & sBoundary) != sBoundary) 404 { 405 areacode |= (sAxis0 & sAxisX) | (sAxis 406 } 407 } 408 return areacode; 409 } 410 else 411 { 412 G4Exception("G4TwistTrapFlatSide::GetAreaC 413 "GeomSolids0001", FatalExcepti 414 "Feature NOT implemented !"); 415 } 416 417 return areacode; 418 } 419 420 //============================================ 421 //* SetCorners ------------------------------- 422 423 void G4TwistTrapFlatSide::SetCorners() 424 { 425 // Set Corner points in local coodinate. 426 427 if (fAxis[0] == kXAxis && fAxis[1] == kYAxi 428 { 429 G4double x, y, z; 430 431 // corner of Axis0min and Axis1min 432 x = -fDx1 + fDy * fTAlph ; 433 y = -fDy ; 434 z = 0 ; 435 SetCorner(sC0Min1Min, x, y, z); 436 437 // corner of Axis0max and Axis1min 438 x = fDx1 + fDy * fTAlph ; 439 y = -fDy ; 440 z = 0 ; 441 SetCorner(sC0Max1Min, x, y, z); 442 443 // corner of Axis0max and Axis1max 444 x = fDx2 - fDy * fTAlph ; 445 y = fDy ; 446 z = 0 ; 447 SetCorner(sC0Max1Max, x, y, z); 448 449 // corner of Axis0min and Axis1max 450 x = -fDx2 - fDy * fTAlph ; 451 y = fDy ; 452 z = 0 ; 453 SetCorner(sC0Min1Max, x, y, z); 454 455 } 456 else 457 { 458 std::ostringstream message; 459 message << "Feature NOT implemented !" << 460 << " fAxis[0] = " << fAxis 461 << " fAxis[1] = " << fAxis 462 G4Exception("G4TwistTrapFlatSide::SetCorn 463 "GeomSolids0001", FatalExcept 464 } 465 } 466 467 //============================================ 468 //* SetBoundaries() -------------------------- 469 470 void G4TwistTrapFlatSide::SetBoundaries() 471 { 472 // Set direction-unit vector of phi-boundar 473 // Don't call the function twice. 474 475 G4ThreeVector direction ; 476 477 if (fAxis[0] == kXAxis && fAxis[1] == kYAxis 478 { 479 // sAxis0 & sAxisMin 480 direction = - ( GetCorner(sC0Min1Max) - Ge 481 direction = direction.unit(); 482 SetBoundary(sAxis0 & (sAxisX | sAxisMin), 483 GetCorner(sC0Min1Max), sAxisY) 484 485 // sAxis0 & sAxisMax 486 direction = GetCorner(sC0Max1Max) - GetCor 487 direction = direction.unit(); 488 SetBoundary(sAxis0 & (sAxisX | sAxisMax), 489 GetCorner(sC0Max1Min), sAxisY) 490 491 // sAxis1 & sAxisMin 492 direction = GetCorner(sC0Max1Min) - GetCor 493 direction = direction.unit(); 494 SetBoundary(sAxis1 & (sAxisY | sAxisMin), 495 GetCorner(sC0Min1Min), sAxisX) 496 497 // sAxis1 & sAxisMax 498 direction = - ( GetCorner(sC0Max1Max) - Ge 499 direction = direction.unit(); 500 SetBoundary(sAxis1 & (sAxisY | sAxisMax), 501 GetCorner(sC0Max1Max), sAxisX) 502 503 } 504 else 505 { 506 std::ostringstream message; 507 message << "Feature NOT implemented !" << 508 << " fAxis[0] = " << fAxis[ 509 << " fAxis[1] = " << fAxis[ 510 G4Exception("G4TwistTrapFlatSide::SetCorne 511 "GeomSolids0001", FatalExcepti 512 } 513 } 514 515 //============================================ 516 //* GetFacets() ------------------------------ 517 518 void G4TwistTrapFlatSide::GetFacets( G4int k, 519 G4int fac 520 { 521 G4double x,y ; // the two parameters 522 G4ThreeVector p ; // a point on the surface 523 524 G4int nnode ; 525 G4int nface ; 526 527 G4double xmin,xmax ; 528 529 // calculate the (n-1)*(k-1) vertices 530 531 for ( G4int i = 0 ; i<n ; ++i ) 532 { 533 y = -fDy + i*(2*fDy)/(n-1) ; 534 535 for ( G4int j = 0 ; j<k ; ++j ) 536 { 537 xmin = GetBoundaryMin(y) ; 538 xmax = GetBoundaryMax(y) ; 539 x = xmin + j*(xmax-xmin)/(k-1) ; 540 541 nnode = GetNode(i,j,k,n,iside) ; 542 p = SurfacePoint(x,y,true) ; // surface 543 544 xyz[nnode][0] = p.x() ; 545 xyz[nnode][1] = p.y() ; 546 xyz[nnode][2] = p.z() ; 547 548 if ( i<n-1 && j<k-1 ) 549 { 550 nface = GetFace(i,j,k,n,iside) ; 551 552 if (fHandedness < 0) // lower side 553 { 554 faces[nface][0] = GetEdgeVisibility( 555 * ( GetNode(i ,j , 556 faces[nface][1] = GetEdgeVisibility( 557 * ( GetNode(i+1,j , 558 faces[nface][2] = GetEdgeVisibility( 559 * ( GetNode(i+1,j+1, 560 faces[nface][3] = GetEdgeVisibility( 561 * ( GetNode(i ,j+1, 562 } 563 else // upper side 564 { 565 faces[nface][0] = GetEdgeVisibility( 566 * ( GetNode(i ,j , 567 faces[nface][1] = GetEdgeVisibility( 568 * ( GetNode(i ,j+1, 569 faces[nface][2] = GetEdgeVisibility( 570 * ( GetNode(i+1,j+1, 571 faces[nface][3] = GetEdgeVisibility( 572 * ( GetNode(i+1,j , 573 } 574 } 575 } 576 } 577 } 578