00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044 #include <iomanip>
00045
00046 #include "G4VTwistSurface.hh"
00047 #include "G4GeometryTolerance.hh"
00048
00049 const G4int G4VTwistSurface::sOutside = 0x00000000;
00050 const G4int G4VTwistSurface::sInside = 0x10000000;
00051 const G4int G4VTwistSurface::sBoundary = 0x20000000;
00052 const G4int G4VTwistSurface::sCorner = 0x40000000;
00053 const G4int G4VTwistSurface::sC0Min1Min = 0x40000101;
00054 const G4int G4VTwistSurface::sC0Max1Min = 0x40000201;
00055 const G4int G4VTwistSurface::sC0Max1Max = 0x40000202;
00056 const G4int G4VTwistSurface::sC0Min1Max = 0x40000102;
00057 const G4int G4VTwistSurface::sAxisMin = 0x00000101;
00058 const G4int G4VTwistSurface::sAxisMax = 0x00000202;
00059 const G4int G4VTwistSurface::sAxisX = 0x00000404;
00060 const G4int G4VTwistSurface::sAxisY = 0x00000808;
00061 const G4int G4VTwistSurface::sAxisZ = 0x00000C0C;
00062 const G4int G4VTwistSurface::sAxisRho = 0x00001010;
00063 const G4int G4VTwistSurface::sAxisPhi = 0x00001414;
00064
00065
00066 const G4int G4VTwistSurface::sAxis0 = 0x0000FF00;
00067 const G4int G4VTwistSurface::sAxis1 = 0x000000FF;
00068 const G4int G4VTwistSurface::sSizeMask = 0x00000303;
00069 const G4int G4VTwistSurface::sAxisMask = 0x0000FCFC;
00070 const G4int G4VTwistSurface::sAreaMask = 0XF0000000;
00071
00072
00073
00074
00075 G4VTwistSurface::G4VTwistSurface(const G4String &name)
00076 : fIsValidNorm(false), fName(name)
00077 {
00078
00079 fAxis[0] = kUndefined;
00080 fAxis[1] = kUndefined;
00081 fAxisMin[0] = kInfinity;
00082 fAxisMin[1] = kInfinity;
00083 fAxisMax[0] = kInfinity;
00084 fAxisMax[1] = kInfinity;
00085 fHandedness = 1;
00086
00087 for (G4int i=0; i<4; i++)
00088 {
00089 fCorners[i].set(kInfinity, kInfinity, kInfinity);
00090 fNeighbours[i] = 0;
00091 }
00092
00093 fCurrentNormal.p.set(kInfinity, kInfinity, kInfinity);
00094
00095 fAmIOnLeftSide.me.set(kInfinity, kInfinity, kInfinity);
00096 fAmIOnLeftSide.vec.set(kInfinity, kInfinity, kInfinity);
00097 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
00098 }
00099
00100 G4VTwistSurface::G4VTwistSurface(const G4String &name,
00101 const G4RotationMatrix &rot,
00102 const G4ThreeVector &tlate,
00103 G4int handedness,
00104 const EAxis axis0 ,
00105 const EAxis axis1 ,
00106 G4double axis0min,
00107 G4double axis1min,
00108 G4double axis0max,
00109 G4double axis1max )
00110 : fIsValidNorm(false), fName(name)
00111 {
00112 fAxis[0] = axis0;
00113 fAxis[1] = axis1;
00114 fAxisMin[0] = axis0min;
00115 fAxisMin[1] = axis1min;
00116 fAxisMax[0] = axis0max;
00117 fAxisMax[1] = axis1max;
00118 fHandedness = handedness;
00119 fRot = rot;
00120 fTrans = tlate;
00121
00122 for (G4int i=0; i<4; i++)
00123 {
00124 fCorners[i].set(kInfinity, kInfinity, kInfinity);
00125 fNeighbours[i] = 0;
00126 }
00127
00128 fCurrentNormal.p.set(kInfinity, kInfinity, kInfinity);
00129
00130 fAmIOnLeftSide.me.set(kInfinity, kInfinity, kInfinity);
00131 fAmIOnLeftSide.vec.set(kInfinity, kInfinity, kInfinity);
00132 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
00133 }
00134
00135
00136
00137
00138 G4VTwistSurface::G4VTwistSurface( __void__& )
00139 : fHandedness(0), fIsValidNorm(false), kCarTolerance(0.),
00140 fName("")
00141 {
00142 fAxis[0] = fAxis[1] = kXAxis;
00143 fAxisMin[0] = fAxisMin[1] = 0.;
00144 fAxisMax[0] = fAxisMax[1] = 0.;
00145 fNeighbours[0] = fNeighbours[1] = fNeighbours[2] = fNeighbours[3] = 0;
00146 }
00147
00148
00149
00150
00151 G4VTwistSurface::~G4VTwistSurface()
00152 {
00153 }
00154
00155
00156
00157
00158 G4int G4VTwistSurface::AmIOnLeftSide(const G4ThreeVector &me,
00159 const G4ThreeVector &vec,
00160 G4bool withtol)
00161 {
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171 static const G4double kAngTolerance
00172 = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
00173
00174 G4RotationMatrix unitrot;
00175 static const G4RotationMatrix rottol = unitrot.rotateZ(0.5*kAngTolerance);
00176 static const G4RotationMatrix invrottol = unitrot.rotateZ(-1.*kAngTolerance);
00177
00178 if (fAmIOnLeftSide.me == me
00179 && fAmIOnLeftSide.vec == vec
00180 && fAmIOnLeftSide.withTol == withtol) {
00181 return fAmIOnLeftSide.amIOnLeftSide;
00182 }
00183
00184 fAmIOnLeftSide.me = me;
00185 fAmIOnLeftSide.vec = vec;
00186 fAmIOnLeftSide.withTol = withtol;
00187
00188 G4ThreeVector met = (G4ThreeVector(me.x(), me.y(), 0.)).unit();
00189 G4ThreeVector vect = (G4ThreeVector(vec.x(), vec.y(), 0.)).unit();
00190
00191 G4ThreeVector ivect = invrottol * vect;
00192 G4ThreeVector rvect = rottol * vect;
00193
00194 G4double metcrossvect = met.x() * vect.y() - met.y() * vect.x();
00195
00196 if (withtol) {
00197 if (met.x() * ivect.y() - met.y() * ivect.x() > 0 &&
00198 metcrossvect >= 0) {
00199 fAmIOnLeftSide.amIOnLeftSide = 1;
00200 } else if (met.x() * rvect.y() - met.y() * rvect.x() < 0 &&
00201 metcrossvect <= 0) {
00202 fAmIOnLeftSide.amIOnLeftSide = -1;
00203 } else {
00204 fAmIOnLeftSide.amIOnLeftSide = 0;
00205 }
00206 } else {
00207 if (metcrossvect > 0) {
00208 fAmIOnLeftSide.amIOnLeftSide = 1;
00209 } else if (metcrossvect < 0 ) {
00210 fAmIOnLeftSide.amIOnLeftSide = -1;
00211 } else {
00212 fAmIOnLeftSide.amIOnLeftSide = 0;
00213 }
00214 }
00215
00216 #ifdef G4TWISTDEBUG
00217 G4cout << " === G4VTwistSurface::AmIOnLeftSide() =============="
00218 << G4endl;
00219 G4cout << " Name , returncode : " << fName << " "
00220 << fAmIOnLeftSide.amIOnLeftSide << G4endl;
00221 G4cout << " me, vec : " << std::setprecision(14) << me
00222 << " " << vec << G4endl;
00223 G4cout << " met, vect : " << met << " " << vect << G4endl;
00224 G4cout << " ivec, rvec : " << ivect << " " << rvect << G4endl;
00225 G4cout << " met x vect : " << metcrossvect << G4endl;
00226 G4cout << " met x ivec : " << met.cross(ivect) << G4endl;
00227 G4cout << " met x rvec : " << met.cross(rvect) << G4endl;
00228 G4cout << " =============================================="
00229 << G4endl;
00230 #endif
00231
00232 return fAmIOnLeftSide.amIOnLeftSide;
00233 }
00234
00235
00236
00237
00238 G4double G4VTwistSurface::DistanceToBoundary(G4int areacode,
00239 G4ThreeVector &xx,
00240 const G4ThreeVector &p)
00241 {
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251 G4ThreeVector d;
00252 G4ThreeVector x0;
00253 G4double dist = kInfinity;
00254 G4int boundarytype;
00255
00256 if (IsAxis0(areacode) && IsAxis1(areacode)) {
00257 std::ostringstream message;
00258 message << "Point is in the corner area." << G4endl
00259 << " Point is in the corner area. This function returns"
00260 << G4endl
00261 << " a direction vector of a boundary line." << G4endl
00262 << " areacode = " << areacode;
00263 G4Exception("G4VTwistSurface::DistanceToBoundary()", "GeomSolids0003",
00264 FatalException, message);
00265 } else if (IsAxis0(areacode) || IsAxis1(areacode)) {
00266 GetBoundaryParameters(areacode, d, x0, boundarytype);
00267 if (boundarytype == sAxisPhi) {
00268 G4double t = x0.getRho() / p.getRho();
00269 xx.set(t*p.x(), t*p.y(), x0.z());
00270 dist = (xx - p).mag();
00271 } else {
00272
00273
00274 dist = DistanceToLine(p, x0, d, xx);
00275 }
00276 } else {
00277 std::ostringstream message;
00278 message << "Bad areacode of boundary." << G4endl
00279 << " areacode = " << areacode;
00280 G4Exception("G4VTwistSurface::DistanceToBoundary()", "GeomSolids0003",
00281 FatalException, message);
00282 }
00283 return dist;
00284 }
00285
00286
00287
00288
00289 G4double G4VTwistSurface::DistanceToIn(const G4ThreeVector &gp,
00290 const G4ThreeVector &gv,
00291 G4ThreeVector &gxxbest)
00292 {
00293 #ifdef G4TWISTDEBUG
00294 G4cout << " ~~~~~ G4VTwistSurface::DistanceToIn(p,v) - Start ~~~~~" << G4endl;
00295 G4cout << " Name : " << fName << G4endl;
00296 G4cout << " gp : " << gp << G4endl;
00297 G4cout << " gv : " << gv << G4endl;
00298 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
00299 #endif
00300
00301 G4ThreeVector gxx[G4VSURFACENXX];
00302 G4double distance[G4VSURFACENXX] ;
00303 G4int areacode[G4VSURFACENXX] ;
00304 G4bool isvalid[G4VSURFACENXX] ;
00305
00306 for (G4int i = 0 ; i<G4VSURFACENXX ; i++ ) {
00307 distance[i] = kInfinity ;
00308 areacode[i] = sOutside ;
00309 isvalid[i] = false ;
00310 }
00311
00312 G4double bestdistance = kInfinity;
00313 #ifdef G4TWISTDEBUG
00314 G4int besti = -1;
00315 #endif
00316 G4ThreeVector bestgxx(kInfinity, kInfinity, kInfinity);
00317
00318 G4int nxx = DistanceToSurface(gp, gv, gxx, distance, areacode,
00319 isvalid, kValidateWithTol);
00320
00321 for (G4int i=0; i< nxx; i++) {
00322
00323
00324
00325
00326
00327 if (!isvalid[i]) {
00328
00329 continue;
00330 }
00331
00332 G4ThreeVector normal = GetNormal(gxx[i], true);
00333
00334 if ((normal * gv) >= 0) {
00335
00336 #ifdef G4TWISTDEBUG
00337 G4cout << " G4VTwistSurface::DistanceToIn(p,v): "
00338 << "particle goes outword the surface." << G4endl;
00339 #endif
00340 continue;
00341 }
00342
00343
00344
00345
00346
00347 if (IsInside(areacode[i])) {
00348 if (distance[i] < bestdistance) {
00349 bestdistance = distance[i];
00350 bestgxx = gxx[i];
00351 #ifdef G4TWISTDEBUG
00352 besti = i;
00353 G4cout << " G4VTwistSurface::DistanceToIn(p,v): "
00354 << " areacode sInside name, distance = "
00355 << fName << " "<< bestdistance << G4endl;
00356 #endif
00357 }
00358
00359
00360
00361
00362
00363 } else {
00364
00365 G4VTwistSurface *neighbours[2];
00366 G4bool isaccepted[2] = {false, false};
00367 G4int nneighbours = GetNeighbours(areacode[i], neighbours);
00368
00369 for (G4int j=0; j< nneighbours; j++) {
00370
00371
00372
00373 G4ThreeVector tmpgxx[G4VSURFACENXX];
00374 G4double tmpdist[G4VSURFACENXX] ;
00375 G4int tmpareacode[G4VSURFACENXX] ;
00376 G4bool tmpisvalid[G4VSURFACENXX] ;
00377
00378 for (G4int l = 0 ; l<G4VSURFACENXX ; l++ ) {
00379 tmpdist[l] = kInfinity ;
00380 tmpareacode[l] = sOutside ;
00381 tmpisvalid[l] = false ;
00382 }
00383
00384 G4int tmpnxx = neighbours[j]->DistanceToSurface(
00385 gp, gv, tmpgxx, tmpdist,
00386 tmpareacode, tmpisvalid,
00387 kValidateWithTol);
00388 G4ThreeVector neighbournormal;
00389
00390 for (G4int k=0; k< tmpnxx; k++) {
00391
00392
00393
00394
00395
00396
00397
00398 if (IsInside(tmpareacode[k])) {
00399
00400 #ifdef G4TWISTDEBUG
00401 G4cout << " G4VTwistSurface:DistanceToIn(p,v): "
00402 << " intersection "<< tmpgxx[k] << G4endl
00403 << " is inside of neighbour surface of " << fName
00404 << " . returning kInfinity." << G4endl;
00405 G4cout << "~~~~~ G4VTwistSurface::DistanceToIn(p,v) - return ~~~~"
00406 << G4endl;
00407 G4cout << " No intersections " << G4endl;
00408 G4cout << " Name : " << fName << G4endl;
00409 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"
00410 << G4endl;
00411 #endif
00412 if (tmpisvalid[k]) return kInfinity;
00413 continue;
00414
00415
00416
00417
00418
00419
00420 } else if (IsSameBoundary(this,areacode[i],
00421 neighbours[j], tmpareacode[k])) {
00422
00423
00424 neighbournormal = neighbours[j]->GetNormal(tmpgxx[k], true);
00425 if (neighbournormal * gv < 0) isaccepted[j] = true;
00426 }
00427 }
00428
00429
00430
00431
00432 if (nneighbours == 1) isaccepted[1] = true;
00433
00434 }
00435
00436
00437
00438 if (isaccepted[0] == true && isaccepted[1] == true) {
00439 if (distance[i] < bestdistance) {
00440 bestdistance = distance[i];
00441 gxxbest = gxx[i];
00442 #ifdef G4TWISTDEBUG
00443 besti = i;
00444 G4cout << " G4VTwistSurface::DistanceToIn(p,v): "
00445 << " areacode sBoundary & sBoundary distance = "
00446 << fName << " " << distance[i] << G4endl;
00447 #endif
00448 }
00449 }
00450
00451 }
00452 }
00453
00454 gxxbest = bestgxx;
00455
00456 #ifdef G4TWISTDEBUG
00457 if (besti < 0) {
00458 G4cout << "~~~~~ G4VTwistSurface::DistanceToIn(p,v) - return ~~~~" << G4endl;
00459 G4cout << " No intersections " << G4endl;
00460 G4cout << " Name : " << fName << G4endl;
00461 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
00462 } else {
00463 G4cout << "~~~~~ G4VTwistSurface::DistanceToIn(p,v) : return ~~~~" << G4endl;
00464 G4cout << " Name, i : " << fName << " , " << besti << G4endl;
00465 G4cout << " gxx[i] : " << gxxbest << G4endl;
00466 G4cout << " bestdist : " << bestdistance << G4endl;
00467 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
00468 }
00469
00470 #endif
00471
00472 return bestdistance;
00473 }
00474
00475
00476
00477
00478 G4double G4VTwistSurface::DistanceToOut(const G4ThreeVector &gp,
00479 const G4ThreeVector &gv,
00480 G4ThreeVector &gxxbest)
00481 {
00482 #ifdef G4TWISTDEBUG
00483 G4cout << "~~~~~ G4VTwistSurface::DistanceToOut(p,v) - Start ~~~~" << G4endl;
00484 G4cout << " Name : " << fName << G4endl;
00485 G4cout << " gp : " << gp << G4endl;
00486 G4cout << " gv : " << gv << G4endl;
00487 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
00488 #endif
00489
00490 G4ThreeVector gxx[G4VSURFACENXX];
00491 G4double distance[G4VSURFACENXX] ;
00492 G4int areacode[G4VSURFACENXX] ;
00493 G4bool isvalid[G4VSURFACENXX] ;
00494 G4int i;
00495
00496 for ( i = 0 ; i<G4VSURFACENXX ; i++ )
00497 {
00498 distance[i] = kInfinity ;
00499 areacode[i] = sOutside ;
00500 isvalid[i] = false ;
00501 }
00502
00503 G4int nxx;
00504 G4double bestdistance = kInfinity;
00505
00506 nxx = DistanceToSurface(gp, gv, gxx, distance, areacode,
00507 isvalid, kValidateWithTol);
00508
00509 for (i=0; i<nxx; i++) {
00510 if (!(isvalid[i])) {
00511 continue;
00512 }
00513
00514 G4ThreeVector normal = GetNormal(gxx[i], true);
00515 if (normal * gv <= 0) {
00516
00517 #ifdef G4TWISTDEBUG
00518 G4cout << " G4VTwistSurface::DistanceToOut(p,v): normal*gv < 0, normal "
00519 << fName << " " << normal
00520 << G4endl;
00521 #endif
00522 } else {
00523
00524 if (distance[i] < bestdistance) {
00525 bestdistance = distance[i];
00526 gxxbest = gxx[i];
00527 }
00528 }
00529 }
00530
00531 #ifdef G4TWISTDEBUG
00532 if (besti < 0) {
00533 G4cout << "~~~~~ G4VTwistSurface::DistanceToOut(p,v) - return ~~~" << G4endl;
00534 G4cout << " No intersections " << G4endl;
00535 G4cout << " Name : " << fName << G4endl;
00536 G4cout << " bestdist : " << bestdistance << G4endl;
00537 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
00538 } else {
00539 G4cout << "~~~~~ G4VTwistSurface::DistanceToOut(p,v) : return ~~~" << G4endl;
00540 G4cout << " Name, i : " << fName << " , " << i << G4endl;
00541 G4cout << " gxx[i] : " << gxxbest << G4endl;
00542 G4cout << " bestdist : " << bestdistance << G4endl;
00543 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
00544 }
00545 #endif
00546
00547 return bestdistance;
00548 }
00549
00550
00551
00552
00553 G4double G4VTwistSurface::DistanceTo(const G4ThreeVector &gp,
00554 G4ThreeVector &gxxbest)
00555 {
00556 #ifdef G4TWISTDEBUG
00557 G4cout << "~~~~~ G4VTwistSurface::DistanceTo(p) - Start ~~~~~~~~~" << G4endl;
00558 G4cout << " Name : " << fName << G4endl;
00559 G4cout << " gp : " << gp << G4endl;
00560 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
00561 #endif
00562
00563
00564 G4ThreeVector gxx[G4VSURFACENXX];
00565 G4double distance[G4VSURFACENXX] ;
00566 G4int areacode[G4VSURFACENXX] ;
00567
00568 for (G4int i = 0 ; i<G4VSURFACENXX ; i++ ) {
00569 distance[i] = kInfinity ;
00570 areacode[i] = sOutside ;
00571 }
00572
00573
00574 DistanceToSurface(gp, gxx, distance, areacode);
00575 gxxbest = gxx[0];
00576
00577 #ifdef G4TWISTDEBUG
00578 G4cout << "~~~~~ G4VTwistSurface::DistanceTo(p) - return ~~~~~~~~" << G4endl;
00579 G4cout << " Name : " << fName << G4endl;
00580 G4cout << " gxx : " << gxxbest << G4endl;
00581 G4cout << " bestdist : " << distance[0] << G4endl;
00582 G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
00583 #endif
00584
00585 return distance[0];
00586 }
00587
00588
00589
00590
00591 G4bool G4VTwistSurface::IsSameBoundary(G4VTwistSurface *surface1, G4int areacode1,
00592 G4VTwistSurface *surface2, G4int areacode2 ) const
00593 {
00594
00595
00596
00597
00598
00599
00600 G4bool testbitmode = true;
00601 G4bool iscorner[2] = {IsCorner(areacode1, testbitmode),
00602 IsCorner(areacode2, testbitmode)};
00603
00604 if (iscorner[0] && iscorner[1]) {
00605
00606 G4ThreeVector corner1 =
00607 surface1->ComputeGlobalPoint(surface1->GetCorner(areacode1));
00608 G4ThreeVector corner2 =
00609 surface2->ComputeGlobalPoint(surface2->GetCorner(areacode2));
00610
00611 if ((corner1 - corner2).mag() < kCarTolerance) {
00612 return true;
00613 } else {
00614 return false;
00615 }
00616
00617 } else if ((IsBoundary(areacode1, testbitmode) && (!iscorner[0])) &&
00618 (IsBoundary(areacode2, testbitmode) && (!iscorner[1]))) {
00619
00620 G4ThreeVector d1, d2, ld1, ld2;
00621 G4ThreeVector x01, x02, lx01, lx02;
00622 G4int type1, type2;
00623 surface1->GetBoundaryParameters(areacode1, ld1, lx01, type1);
00624 surface2->GetBoundaryParameters(areacode2, ld2, lx02, type2);
00625
00626 x01 = surface1->ComputeGlobalPoint(lx01);
00627 x02 = surface2->ComputeGlobalPoint(lx02);
00628 d1 = surface1->ComputeGlobalDirection(ld1);
00629 d2 = surface2->ComputeGlobalDirection(ld2);
00630
00631 if ((x01 - x02).mag() < kCarTolerance &&
00632 (d1 - d2).mag() < kCarTolerance) {
00633 return true;
00634 } else {
00635 return false;
00636 }
00637
00638 } else {
00639 return false;
00640 }
00641 }
00642
00643
00644
00645
00646 void G4VTwistSurface::GetBoundaryParameters(const G4int &areacode,
00647 G4ThreeVector &d,
00648 G4ThreeVector &x0,
00649 G4int &boundarytype) const
00650 {
00651
00652
00653
00654
00655 G4int i;
00656 for (i=0; i<4; i++) {
00657 if (fBoundaries[i].GetBoundaryParameters(areacode, d, x0,
00658 boundarytype)) {
00659 return;
00660 }
00661 }
00662
00663 std::ostringstream message;
00664 message << "Not registered boundary." << G4endl
00665 << " Boundary at areacode " << std::hex << areacode
00666 << std::dec << G4endl
00667 << " is not registered.";
00668 G4Exception("G4VTwistSurface::GetBoundaryParameters()", "GeomSolids0002",
00669 FatalException, message);
00670 }
00671
00672
00673
00674
00675 G4ThreeVector G4VTwistSurface::GetBoundaryAtPZ(G4int areacode,
00676 const G4ThreeVector &p) const
00677 {
00678
00679
00680
00681
00682 if (areacode & sAxis0 && areacode & sAxis1) {
00683 std::ostringstream message;
00684 message << "Point is in the corner area." << G4endl
00685 << " This function returns "
00686 << "a direction vector of a boundary line." << G4endl
00687 << " areacode = " << areacode;
00688 G4Exception("G4VTwistSurface::GetBoundaryAtPZ()", "GeomSolids0003",
00689 FatalException, message);
00690 }
00691
00692 G4ThreeVector d;
00693 G4ThreeVector x0;
00694 G4int boundarytype;
00695 G4bool found = false;
00696
00697 for (G4int i=0; i<4; i++) {
00698 if (fBoundaries[i].GetBoundaryParameters(areacode, d, x0,
00699 boundarytype)){
00700 found = true;
00701 continue;
00702 }
00703 }
00704
00705 if (!found) {
00706 std::ostringstream message;
00707 message << "Not registered boundary." << G4endl
00708 << " Boundary at areacode " << areacode << G4endl
00709 << " is not registered.";
00710 G4Exception("G4VTwistSurface::GetBoundaryAtPZ()", "GeomSolids0002",
00711 FatalException, message);
00712 }
00713
00714 if (((boundarytype & sAxisPhi) == sAxisPhi) ||
00715 ((boundarytype & sAxisRho) == sAxisRho)) {
00716 std::ostringstream message;
00717 message << "Not a z-depended line boundary." << G4endl
00718 << " Boundary at areacode " << areacode << G4endl
00719 << " is not a z-depended line.";
00720 G4Exception("G4VTwistSurface::GetBoundaryAtPZ()", "GeomSolids0002",
00721 FatalException, message);
00722 }
00723 return ((p.z() - x0.z()) / d.z()) * d + x0;
00724 }
00725
00726
00727
00728
00729 void G4VTwistSurface::SetCorner(G4int areacode, G4double x, G4double y, G4double z)
00730 {
00731 if ((areacode & sCorner) != sCorner){
00732 std::ostringstream message;
00733 message << "Area code must represents corner." << G4endl
00734 << " areacode " << areacode;
00735 G4Exception("G4VTwistSurface::SetCorner()", "GeomSolids0002",
00736 FatalException, message);
00737 }
00738
00739 if ((areacode & sC0Min1Min) == sC0Min1Min) {
00740 fCorners[0].set(x, y, z);
00741 } else if ((areacode & sC0Max1Min) == sC0Max1Min) {
00742 fCorners[1].set(x, y, z);
00743 } else if ((areacode & sC0Max1Max) == sC0Max1Max) {
00744 fCorners[2].set(x, y, z);
00745 } else if ((areacode & sC0Min1Max) == sC0Min1Max) {
00746 fCorners[3].set(x, y, z);
00747 }
00748 }
00749
00750
00751
00752
00753 void G4VTwistSurface::GetBoundaryAxis(G4int areacode, EAxis axis[]) const
00754 {
00755 if ((areacode & sBoundary) != sBoundary) {
00756 G4Exception("G4VTwistSurface::GetBoundaryAxis()", "GeomSolids0003",
00757 FatalException, "Not located on a boundary!");
00758 }
00759 G4int i;
00760 for (i=0; i<2; i++) {
00761
00762 G4int whichaxis = 0 ;
00763 if (i == 0) {
00764 whichaxis = sAxis0;
00765 } else if (i == 1) {
00766 whichaxis = sAxis1;
00767 }
00768
00769
00770 G4int axiscode = whichaxis & sAxisMask & areacode ;
00771 if (axiscode) {
00772 if (axiscode == (whichaxis & sAxisX)) {
00773 axis[i] = kXAxis;
00774 } else if (axiscode == (whichaxis & sAxisY)) {
00775 axis[i] = kYAxis;
00776 } else if (axiscode == (whichaxis & sAxisZ)) {
00777 axis[i] = kZAxis;
00778 } else if (axiscode == (whichaxis & sAxisRho)) {
00779 axis[i] = kRho;
00780 } else if (axiscode == (whichaxis & sAxisPhi)) {
00781 axis[i] = kPhi;
00782 } else {
00783 std::ostringstream message;
00784 message << "Not supported areacode." << G4endl
00785 << " areacode " << areacode;
00786 G4Exception("G4VTwistSurface::GetBoundaryAxis()", "GeomSolids0001",
00787 FatalException, message);
00788 }
00789 }
00790 }
00791 }
00792
00793
00794
00795
00796 void G4VTwistSurface::GetBoundaryLimit(G4int areacode, G4double limit[]) const
00797 {
00798 if (areacode & sCorner) {
00799 if (areacode & sC0Min1Min) {
00800 limit[0] = fAxisMin[0];
00801 limit[1] = fAxisMin[1];
00802 } else if (areacode & sC0Max1Min) {
00803 limit[0] = fAxisMax[0];
00804 limit[1] = fAxisMin[1];
00805 } else if (areacode & sC0Max1Max) {
00806 limit[0] = fAxisMax[0];
00807 limit[1] = fAxisMax[1];
00808 } else if (areacode & sC0Min1Max) {
00809 limit[0] = fAxisMin[0];
00810 limit[1] = fAxisMax[1];
00811 }
00812 } else if (areacode & sBoundary) {
00813 if (areacode & (sAxis0 | sAxisMin)) {
00814 limit[0] = fAxisMin[0];
00815 } else if (areacode & (sAxis1 | sAxisMin)) {
00816 limit[0] = fAxisMin[1];
00817 } else if (areacode & (sAxis0 | sAxisMax)) {
00818 limit[0] = fAxisMax[0];
00819 } else if (areacode & (sAxis1 | sAxisMax)) {
00820 limit[0] = fAxisMax[1];
00821 }
00822 } else {
00823 std::ostringstream message;
00824 message << "Not located on a boundary!" << G4endl
00825 << " areacode " << areacode;
00826 G4Exception("G4VTwistSurface::GetBoundaryLimit()", "GeomSolids1002",
00827 JustWarning, message);
00828 }
00829 }
00830
00831
00832
00833
00834 void G4VTwistSurface::SetBoundary(const G4int &axiscode,
00835 const G4ThreeVector &direction,
00836 const G4ThreeVector &x0,
00837 const G4int &boundarytype)
00838 {
00839 G4int code = (~sAxisMask) & axiscode;
00840 if ((code == (sAxis0 & sAxisMin)) ||
00841 (code == (sAxis0 & sAxisMax)) ||
00842 (code == (sAxis1 & sAxisMin)) ||
00843 (code == (sAxis1 & sAxisMax))) {
00844
00845 G4int i;
00846 G4bool done = false;
00847 for (i=0; i<4; i++) {
00848 if (fBoundaries[i].IsEmpty()) {
00849 fBoundaries[i].SetFields(axiscode, direction,
00850 x0, boundarytype);
00851 done = true;
00852 break;
00853 }
00854 }
00855
00856 if (!done) {
00857 G4Exception("G4VTwistSurface::SetBoundary()", "GeomSolids0003",
00858 FatalException, "Number of boundary exceeding 4!");
00859 }
00860 } else {
00861 std::ostringstream message;
00862 message << "Invalid axis-code." << G4endl
00863 << " axiscode = "
00864 << std::hex << axiscode << std::dec;
00865 G4Exception("G4VTwistSurface::SetBoundary()", "GeomSolids0003",
00866 FatalException, message);
00867 }
00868 }
00869
00870
00871
00872
00873 G4int G4VTwistSurface::GetFace( G4int i, G4int j, G4int k,
00874 G4int n, G4int iside )
00875 {
00876
00877
00878
00879 if ( iside == 0 ) {
00880 return i * ( k - 1 ) + j ;
00881 }
00882
00883 else if ( iside == 1 ) {
00884 return (k-1)*(k-1) + i*(k-1) + j ;
00885 }
00886
00887 else if ( iside == 2 ) {
00888 return 2*(k-1)*(k-1) + i*(k-1) + j ;
00889 }
00890
00891 else if ( iside == 3 ) {
00892 return 2*(k-1)*(k-1) + (n-1)*(k-1) + i*(k-1) + j ;
00893 }
00894
00895 else if ( iside == 4 ) {
00896 return 2*(k-1)*(k-1) + 2*(n-1)*(k-1) + i*(k-1) + j ;
00897 }
00898
00899 else if ( iside == 5 ) {
00900 return 2*(k-1)*(k-1) + 3*(n-1)*(k-1) + i*(k-1) + j ;
00901 }
00902
00903 else {
00904
00905 std::ostringstream message;
00906 message << "Not correct side number: "
00907 << GetName() << G4endl
00908 << "iside is " << iside << " but should be "
00909 << "0,1,2,3,4 or 5" << ".";
00910 G4Exception("G4TwistSurface::G4GetFace()", "GeomSolids0002",
00911 FatalException, message);
00912
00913
00914 }
00915
00916 return -1 ;
00917 }
00918
00919
00920
00921
00922 G4int G4VTwistSurface::GetNode( G4int i, G4int j, G4int k,
00923 G4int n, G4int iside )
00924 {
00925
00926
00927
00928
00929 if ( iside == 0 ) {
00930
00931
00932 return i * k + j ;
00933 }
00934
00935 if ( iside == 1 ) {
00936
00937
00938 return k*k + i*k + j ;
00939 }
00940
00941 else if ( iside == 2 ) {
00942
00943 if ( i == 0 ) { return j ; }
00944 else if ( i == n-1 ) { return k*k + j ; }
00945 else { return 2*k*k + 4*(i-1)*(k-1) + j ; }
00946 }
00947
00948 else if ( iside == 3 ) {
00949
00950 if ( i == 0 ) { return (j+1)*k - 1 ; }
00951 else if ( i == n-1 ) { return k*k + (j+1)*k - 1 ; }
00952 else { return 2*k*k + 4*(i-1)*(k-1) + (k-1) + j ; }
00953 }
00954 else if ( iside == 4 ) {
00955
00956 if ( i == 0 ) { return k*k - 1 - j ; }
00957 else if ( i == n-1 ) { return 2*k*k - 1 - j ; }
00958 else { return 2*k*k + 4*(i-1)*(k-1) + 2*(k-1) + j ;
00959 }
00960 }
00961 else if ( iside == 5 ) {
00962
00963 if ( i == 0 ) { return k*k - (j+1)*k ; }
00964 else if ( i == n-1) { return 2*k*k - (j+1)*k ; }
00965 else {
00966 if ( j == k-1 ) { return 2*k*k + 4*(i-1)*(k-1) ; }
00967 else { return 2*k*k + 4*(i-1)*(k-1) + 3*(k-1) + j ; }
00968 }
00969 }
00970
00971 else {
00972
00973 std::ostringstream message;
00974 message << "Not correct side number: "
00975 << GetName() << G4endl
00976 << "iside is " << iside << " but should be "
00977 << "0,1,2,3,4 or 5" << ".";
00978 G4Exception("G4TwistSurface::G4GetNode()", "GeomSolids0002",
00979 FatalException, message);
00980 }
00981 return -1 ;
00982 }
00983
00984
00985
00986
00987 G4int G4VTwistSurface::GetEdgeVisibility( G4int i, G4int j, G4int k, G4int n, G4int number, G4int orientation)
00988 {
00989
00990
00991
00992
00993
00994
00995
00996
00997
00998
00999
01000
01001
01002
01003
01004
01005
01006
01007
01008
01009
01010
01011
01012
01013
01014 if ( ( i>0 && i<n-2 ) && ( j>0 && j<k-2 ) ) {
01015 return -1 ;
01016 }
01017
01018
01019
01020 if ( orientation < 0 ) { number = ( 3 - number ) ; }
01021
01022
01023 if ( ( j>=1 && j<=k-3 ) ) {
01024
01025 if ( i == 0 ) {
01026 return ( number == 3 ) ? 1 : -1 ;
01027 }
01028
01029 else if ( i == n-2 ) {
01030 return ( number == 1 ) ? 1 : -1 ;
01031 }
01032
01033 else {
01034 std::ostringstream message;
01035 message << "Not correct face number: " << GetName() << " !";
01036 G4Exception("G4TwistSurface::G4GetEdgeVisibility()",
01037 "GeomSolids0003", FatalException, message);
01038 }
01039 }
01040
01041 if ( ( i>=1 && i<=n-3 ) ) {
01042
01043 if ( j == 0 ) {
01044 return ( number == 0 ) ? 1 : -1 ;
01045 }
01046
01047 else if ( j == k-2 ) {
01048 return ( number == 2 ) ? 1 : -1 ;
01049 }
01050
01051 else {
01052 std::ostringstream message;
01053 message << "Not correct face number: " << GetName() << " !";
01054 G4Exception("G4TwistSurface::G4GetEdgeVisibility()",
01055 "GeomSolids0003", FatalException, message);
01056 }
01057 }
01058
01059
01060 if ( i == 0 && j == 0 ) {
01061 return ( number == 0 || number == 3 ) ? 1 : -1 ;
01062 }
01063 else if ( i == 0 && j == k-2 ) {
01064 return ( number == 2 || number == 3 ) ? 1 : -1 ;
01065 }
01066 else if ( i == n-2 && j == k-2 ) {
01067 return ( number == 1 || number == 2 ) ? 1 : -1 ;
01068 }
01069 else if ( i == n-2 && j == 0 ) {
01070 return ( number == 0 || number == 1 ) ? 1 : -1 ;
01071 }
01072 else {
01073 std::ostringstream message;
01074 message << "Not correct face number: " << GetName() << " !";
01075 G4Exception("G4TwistSurface::G4GetEdgeVisibility()",
01076 "GeomSolids0003", FatalException, message);
01077 }
01078
01079 std::ostringstream message;
01080 message << "Not correct face number: " << GetName() << " !";
01081 G4Exception("G4TwistSurface::G4GetEdgeVisibility()", "GeomSolids0003",
01082 FatalException, message);
01083
01084 return 0 ;
01085 }
01086
01087
01088
01089
01090
01091 void G4VTwistSurface::DebugPrint() const
01092 {
01093 G4ThreeVector A = fRot * GetCorner(sC0Min1Min) + fTrans;
01094 G4ThreeVector B = fRot * GetCorner(sC0Max1Min) + fTrans;
01095 G4ThreeVector C = fRot * GetCorner(sC0Max1Max) + fTrans;
01096 G4ThreeVector D = fRot * GetCorner(sC0Min1Max) + fTrans;
01097
01098 G4cout << "/* G4VTwistSurface::DebugPrint():-------------------------------"
01099 << G4endl;
01100 G4cout << "/* Name = " << fName << G4endl;
01101 G4cout << "/* Axis = " << std::hex << fAxis[0] << " "
01102 << std::hex << fAxis[1]
01103 << " (0,1,2,3,5 = kXAxis,kYAxis,kZAxis,kRho,kPhi)"
01104 << std::dec << G4endl;
01105 G4cout << "/* BoundaryLimit(in local) fAxis0(min, max) = ("<<fAxisMin[0]
01106 << ", " << fAxisMax[0] << ")" << G4endl;
01107 G4cout << "/* BoundaryLimit(in local) fAxis1(min, max) = ("<<fAxisMin[1]
01108 << ", " << fAxisMax[1] << ")" << G4endl;
01109 G4cout << "/* Cornar point sC0Min1Min = " << A << G4endl;
01110 G4cout << "/* Cornar point sC0Max1Min = " << B << G4endl;
01111 G4cout << "/* Cornar point sC0Max1Max = " << C << G4endl;
01112 G4cout << "/* Cornar point sC0Min1Max = " << D << G4endl;
01113 G4cout << "/*---------------------------------------------------------"
01114 << G4endl;
01115 }
01116
01117
01118
01119
01120
01121
01122
01123
01124 G4VTwistSurface::CurrentStatus::CurrentStatus()
01125 {
01126 for (size_t i=0; i<G4VSURFACENXX; i++)
01127 {
01128 fDistance[i] = kInfinity;
01129 fAreacode[i] = sOutside;
01130 fIsValid[i] = false;
01131 fXX[i].set(kInfinity, kInfinity, kInfinity);
01132 }
01133 fNXX = 0;
01134 fLastp.set(kInfinity, kInfinity, kInfinity);
01135 fLastv.set(kInfinity, kInfinity, kInfinity);
01136 fLastValidate = kUninitialized;
01137 fDone = false;
01138 }
01139
01140
01141
01142
01143 G4VTwistSurface::CurrentStatus::~CurrentStatus()
01144 {
01145 }
01146
01147
01148
01149
01150 void
01151 G4VTwistSurface::CurrentStatus::SetCurrentStatus(G4int i,
01152 G4ThreeVector &xx,
01153 G4double &dist,
01154 G4int &areacode,
01155 G4bool &isvalid,
01156 G4int nxx,
01157 EValidate validate,
01158 const G4ThreeVector *p,
01159 const G4ThreeVector *v)
01160 {
01161 fDistance[i] = dist;
01162 fAreacode[i] = areacode;
01163 fIsValid[i] = isvalid;
01164 fXX[i] = xx;
01165 fNXX = nxx;
01166 fLastValidate = validate;
01167 if (p)
01168 {
01169 fLastp = *p;
01170 }
01171 else
01172 {
01173 G4Exception("G4VTwistSurface::CurrentStatus::SetCurrentStatus()",
01174 "GeomSolids0003", FatalException, "SetCurrentStatus: p = 0!");
01175 }
01176 if (v)
01177 {
01178 fLastv = *v;
01179 }
01180 else
01181 {
01182 fLastv.set(kInfinity, kInfinity, kInfinity);
01183 }
01184 fDone = true;
01185 }
01186
01187
01188
01189
01190 void
01191 G4VTwistSurface::CurrentStatus::ResetfDone(EValidate validate,
01192 const G4ThreeVector *p,
01193 const G4ThreeVector *v)
01194
01195 {
01196 if (validate == fLastValidate && p && *p == fLastp)
01197 {
01198 if (!v || (v && *v == fLastv)) return;
01199 }
01200 G4ThreeVector xx(kInfinity, kInfinity, kInfinity);
01201 for (size_t i=0; i<G4VSURFACENXX; i++)
01202 {
01203 fDistance[i] = kInfinity;
01204 fAreacode[i] = sOutside;
01205 fIsValid[i] = false;
01206 fXX[i] = xx;
01207 }
01208 fNXX = 0;
01209 fLastp.set(kInfinity, kInfinity, kInfinity);
01210 fLastv.set(kInfinity, kInfinity, kInfinity);
01211 fLastValidate = kUninitialized;
01212 fDone = false;
01213 }
01214
01215
01216
01217
01218 void
01219 G4VTwistSurface::CurrentStatus::DebugPrint() const
01220 {
01221 G4cout << "CurrentStatus::Dist0,1= " << fDistance[0]
01222 << " " << fDistance[1] << " areacode = " << fAreacode[0]
01223 << " " << fAreacode[1] << G4endl;
01224 }
01225
01226
01227
01228
01229
01230
01231
01232
01233 G4VTwistSurface::Boundary::Boundary()
01234 : fBoundaryAcode(-1), fBoundaryType(0)
01235 {
01236 }
01237
01238
01239
01240
01241 G4VTwistSurface::Boundary::~Boundary()
01242 {
01243 }
01244
01245
01246
01247
01248 void
01249 G4VTwistSurface::Boundary::SetFields(const G4int &areacode,
01250 const G4ThreeVector &d,
01251 const G4ThreeVector &x0,
01252 const G4int &boundarytype)
01253 {
01254 fBoundaryAcode = areacode;
01255 fBoundaryDirection = d;
01256 fBoundaryX0 = x0;
01257 fBoundaryType = boundarytype;
01258 }
01259
01260
01261
01262
01263 G4bool G4VTwistSurface::Boundary::IsEmpty() const
01264 {
01265 if (fBoundaryAcode == -1) return true;
01266 return false;
01267 }
01268
01269
01270
01271
01272 G4bool
01273 G4VTwistSurface::Boundary::GetBoundaryParameters(const G4int &areacode,
01274 G4ThreeVector &d,
01275 G4ThreeVector &x0,
01276 G4int &boundarytype) const
01277 {
01278
01279
01280
01281 if ((areacode & sAxis0) && (areacode & sAxis1))
01282 {
01283 std::ostringstream message;
01284 message << "Located in the corner area." << G4endl
01285 << " This function returns a direction vector of "
01286 << "a boundary line." << G4endl
01287 << " areacode = " << areacode;
01288 G4Exception("G4VTwistSurface::Boundary::GetBoundaryParameters()",
01289 "GeomSolids0003", FatalException, message);
01290 }
01291 if ((areacode & sSizeMask) != (fBoundaryAcode & sSizeMask))
01292 {
01293 return false;
01294 }
01295 d = fBoundaryDirection;
01296 x0 = fBoundaryX0;
01297 boundarytype = fBoundaryType;
01298 return true;
01299 }