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
00045
00046
00047
00049
00050 #include "G4Trap.hh"
00051 #include "globals.hh"
00052
00053 #include "G4VoxelLimits.hh"
00054 #include "G4AffineTransform.hh"
00055
00056 #include "G4VPVParameterisation.hh"
00057
00058 #include "Randomize.hh"
00059
00060 #include "G4VGraphicsScene.hh"
00061 #include "G4Polyhedron.hh"
00062 #include "G4NURBS.hh"
00063 #include "G4NURBSbox.hh"
00064
00065 using namespace CLHEP;
00066
00068
00069
00070
00071 const G4double kCoplanar_Tolerance = 1E-4 ;
00072
00074
00075
00076
00077 enum Eside {kUndef,ks0,ks1,ks2,ks3,kPZ,kMZ};
00078
00080
00081
00082
00083
00084 G4Trap::G4Trap( const G4String& pName,
00085 G4double pDz,
00086 G4double pTheta, G4double pPhi,
00087 G4double pDy1, G4double pDx1, G4double pDx2,
00088 G4double pAlp1,
00089 G4double pDy2, G4double pDx3, G4double pDx4,
00090 G4double pAlp2)
00091 : G4CSGSolid(pName)
00092 {
00093 if ( pDz <= 0 || pDy1 <= 0 || pDx1 <= 0 ||
00094 pDx2 <= 0 || pDy2 <= 0 || pDx3 <= 0 || pDx4 <= 0 )
00095 {
00096 std::ostringstream message;
00097 message << "Invalid length parameters for Solid: " << GetName() << G4endl
00098 << " X - "
00099 << pDx1 << ", " << pDx2 << ", " << pDx3 << ", " << pDx4 << G4endl
00100 << " Y - " << pDy1 << ", " << pDy2 << G4endl
00101 << " Z - " << pDz;
00102 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
00103 FatalException, message);
00104 }
00105
00106 fDz=pDz;
00107 fTthetaCphi=std::tan(pTheta)*std::cos(pPhi);
00108 fTthetaSphi=std::tan(pTheta)*std::sin(pPhi);
00109
00110 fDy1=pDy1;
00111 fDx1=pDx1;
00112 fDx2=pDx2;
00113 fTalpha1=std::tan(pAlp1);
00114
00115 fDy2=pDy2;
00116 fDx3=pDx3;
00117 fDx4=pDx4;
00118 fTalpha2=std::tan(pAlp2);
00119
00120 MakePlanes();
00121 }
00122
00124
00125
00126
00127
00128
00129 G4Trap::G4Trap( const G4String& pName,
00130 const G4ThreeVector pt[8] )
00131 : G4CSGSolid(pName)
00132 {
00133 G4bool good;
00134
00135
00136
00137
00138 if (!( pt[0].z() < 0
00139 && pt[0].z() == pt[1].z() && pt[0].z() == pt[2].z()
00140 && pt[0].z() == pt[3].z()
00141 && pt[4].z() > 0
00142 && pt[4].z() == pt[5].z() && pt[4].z() == pt[6].z()
00143 && pt[4].z() == pt[7].z()
00144 && std::fabs( pt[0].z() + pt[4].z() ) < kCarTolerance
00145 && pt[0].y() == pt[1].y() && pt[2].y() == pt[3].y()
00146 && pt[4].y() == pt[5].y() && pt[6].y() == pt[7].y()
00147 && std::fabs( pt[0].y() + pt[2].y() + pt[4].y() + pt[6].y() ) < kCarTolerance
00148 && std::fabs( pt[0].x() + pt[1].x() + pt[4].x() + pt[5].x() +
00149 pt[2].x() + pt[3].x() + pt[6].x() + pt[7].x() ) < kCarTolerance ) )
00150 {
00151 std::ostringstream message;
00152 message << "Invalid vertice coordinates for Solid: " << GetName();
00153 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
00154 FatalException, message);
00155 }
00156
00157
00158
00159 good = MakePlane(pt[0],pt[4],pt[5],pt[1],fPlanes[0]);
00160
00161 if (!good)
00162 {
00163 DumpInfo();
00164 G4Exception("G4Trap::G4Trap()", "GeomSolids0002", FatalException,
00165 "Face at ~-Y not planar.");
00166 }
00167
00168
00169
00170 good = MakePlane(pt[2],pt[3],pt[7],pt[6],fPlanes[1]);
00171
00172 if (!good)
00173 {
00174 std::ostringstream message;
00175 message << "Face at ~+Y not planar for Solid: " << GetName();
00176 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
00177 FatalException, message);
00178 }
00179
00180
00181
00182 good = MakePlane(pt[0],pt[2],pt[6],pt[4],fPlanes[2]);
00183
00184 if (!good)
00185 {
00186 std::ostringstream message;
00187 message << "Face at ~-X not planar for Solid: " << GetName();
00188 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
00189 FatalException, message);
00190 }
00191
00192
00193
00194 good = MakePlane(pt[1],pt[5],pt[7],pt[3],fPlanes[3]);
00195 if (!good)
00196 {
00197 std::ostringstream message;
00198 message << "Face at ~+X not planar for Solid: " << GetName();
00199 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
00200 FatalException, message);
00201 }
00202 fDz = (pt[7]).z() ;
00203
00204 fDy1 = ((pt[2]).y()-(pt[1]).y())*0.5;
00205 fDx1 = ((pt[1]).x()-(pt[0]).x())*0.5;
00206 fDx2 = ((pt[3]).x()-(pt[2]).x())*0.5;
00207 fTalpha1 = ((pt[2]).x()+(pt[3]).x()-(pt[1]).x()-(pt[0]).x())*0.25/fDy1;
00208
00209 fDy2 = ((pt[6]).y()-(pt[5]).y())*0.5;
00210 fDx3 = ((pt[5]).x()-(pt[4]).x())*0.5;
00211 fDx4 = ((pt[7]).x()-(pt[6]).x())*0.5;
00212 fTalpha2 = ((pt[6]).x()+(pt[7]).x()-(pt[5]).x()-(pt[4]).x())*0.25/fDy2;
00213
00214 fTthetaCphi = ((pt[4]).x()+fDy2*fTalpha2+fDx3)/fDz;
00215 fTthetaSphi = ((pt[4]).y()+fDy2)/fDz;
00216 }
00217
00219
00220
00221
00222 G4Trap::G4Trap( const G4String& pName,
00223 G4double pZ,
00224 G4double pY,
00225 G4double pX, G4double pLTX )
00226 : G4CSGSolid(pName)
00227 {
00228 G4bool good;
00229
00230 if ( pZ<=0 || pY<=0 || pX<=0 || pLTX<=0 || pLTX>pX )
00231 {
00232 std::ostringstream message;
00233 message << "Invalid length parameters for Solid: " << GetName();
00234 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
00235 FatalException, message);
00236 }
00237
00238 fDz = 0.5*pZ ;
00239 fTthetaCphi = 0 ;
00240 fTthetaSphi = 0 ;
00241
00242 fDy1 = 0.5*pY;
00243 fDx1 = 0.5*pX ;
00244 fDx2 = 0.5*pLTX;
00245 fTalpha1 = 0.5*(pLTX - pX)/pY;
00246
00247 fDy2 = fDy1 ;
00248 fDx3 = fDx1;
00249 fDx4 = fDx2 ;
00250 fTalpha2 = fTalpha1 ;
00251
00252 G4ThreeVector pt[8] ;
00253
00254 pt[0]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
00255 -fDz*fTthetaSphi-fDy1,-fDz);
00256 pt[1]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
00257 -fDz*fTthetaSphi-fDy1,-fDz);
00258 pt[2]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
00259 -fDz*fTthetaSphi+fDy1,-fDz);
00260 pt[3]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
00261 -fDz*fTthetaSphi+fDy1,-fDz);
00262 pt[4]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
00263 +fDz*fTthetaSphi-fDy2,+fDz);
00264 pt[5]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
00265 +fDz*fTthetaSphi-fDy2,+fDz);
00266 pt[6]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
00267 +fDz*fTthetaSphi+fDy2,+fDz);
00268 pt[7]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
00269 +fDz*fTthetaSphi+fDy2,+fDz);
00270
00271
00272
00273 good=MakePlane(pt[0],pt[4],pt[5],pt[1],fPlanes[0]);
00274 if (!good)
00275 {
00276 std::ostringstream message;
00277 message << "Face at ~-Y not planar for Solid: " << GetName();
00278 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
00279 FatalException, message);
00280 }
00281
00282
00283
00284 good=MakePlane(pt[2],pt[3],pt[7],pt[6],fPlanes[1]);
00285 if (!good)
00286 {
00287 std::ostringstream message;
00288 message << "Face at ~+Y not planar for Solid: " << GetName();
00289 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
00290 FatalException, message);
00291 }
00292
00293
00294
00295 good=MakePlane(pt[0],pt[2],pt[6],pt[4],fPlanes[2]);
00296 if (!good)
00297 {
00298 std::ostringstream message;
00299 message << "Face at ~-X not planar for Solid: " << GetName();
00300 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
00301 FatalException, message);
00302 }
00303
00304
00305
00306 good=MakePlane(pt[1],pt[5],pt[7],pt[3],fPlanes[3]);
00307 if (!good)
00308 {
00309 std::ostringstream message;
00310 message << "Face at ~+X not planar for Solid: " << GetName();
00311 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
00312 FatalException, message);
00313 }
00314 }
00315
00317
00318
00319
00320 G4Trap::G4Trap( const G4String& pName,
00321 G4double pDx1, G4double pDx2,
00322 G4double pDy1, G4double pDy2,
00323 G4double pDz )
00324 : G4CSGSolid(pName)
00325 {
00326 G4bool good;
00327
00328 if ( pDz<=0 || pDy1<=0 || pDx1<=0 || pDx2<=0 || pDy2<=0 )
00329 {
00330 std::ostringstream message;
00331 message << "Invalid length parameters for Solid: " << GetName();
00332 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
00333 FatalException, message);
00334 }
00335
00336 fDz = pDz;
00337 fTthetaCphi = 0 ;
00338 fTthetaSphi = 0 ;
00339
00340 fDy1 = pDy1 ;
00341 fDx1 = pDx1 ;
00342 fDx2 = pDx1 ;
00343 fTalpha1 = 0 ;
00344
00345 fDy2 = pDy2 ;
00346 fDx3 = pDx2 ;
00347 fDx4 = pDx2 ;
00348 fTalpha2 = 0 ;
00349
00350 G4ThreeVector pt[8] ;
00351
00352 pt[0]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
00353 -fDz*fTthetaSphi-fDy1,-fDz);
00354 pt[1]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
00355 -fDz*fTthetaSphi-fDy1,-fDz);
00356 pt[2]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
00357 -fDz*fTthetaSphi+fDy1,-fDz);
00358 pt[3]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
00359 -fDz*fTthetaSphi+fDy1,-fDz);
00360 pt[4]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
00361 +fDz*fTthetaSphi-fDy2,+fDz);
00362 pt[5]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
00363 +fDz*fTthetaSphi-fDy2,+fDz);
00364 pt[6]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
00365 +fDz*fTthetaSphi+fDy2,+fDz);
00366 pt[7]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
00367 +fDz*fTthetaSphi+fDy2,+fDz);
00368
00369
00370
00371 good=MakePlane(pt[0],pt[4],pt[5],pt[1],fPlanes[0]);
00372 if (!good)
00373 {
00374 std::ostringstream message;
00375 message << "Face at ~-Y not planar for Solid: " << GetName();
00376 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
00377 FatalException, message);
00378 }
00379
00380
00381
00382 good=MakePlane(pt[2],pt[3],pt[7],pt[6],fPlanes[1]);
00383 if (!good)
00384 {
00385 std::ostringstream message;
00386 message << "Face at ~+Y not planar for Solid: " << GetName();
00387 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
00388 FatalException, message);
00389 }
00390
00391
00392
00393 good=MakePlane(pt[0],pt[2],pt[6],pt[4],fPlanes[2]);
00394 if (!good)
00395 {
00396 std::ostringstream message;
00397 message << "Face at ~-X not planar for Solid: " << GetName();
00398 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
00399 FatalException, message);
00400 }
00401
00402
00403
00404 good=MakePlane(pt[1],pt[5],pt[7],pt[3],fPlanes[3]);
00405 if (!good)
00406 {
00407 std::ostringstream message;
00408 message << "Face at ~+X not planar for Solid: " << GetName();
00409 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
00410 FatalException, message);
00411 }
00412 }
00413
00415
00416
00417
00418 G4Trap::G4Trap( const G4String& pName,
00419 G4double pDx, G4double pDy,
00420 G4double pDz,
00421 G4double pAlpha,
00422 G4double pTheta, G4double pPhi )
00423 : G4CSGSolid(pName)
00424 {
00425 G4bool good;
00426
00427 if ( pDz<=0 || pDy<=0 || pDx<=0 )
00428 {
00429 std::ostringstream message;
00430 message << "Invalid length parameters for Solid: " << GetName();
00431 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
00432 FatalException, message);
00433 }
00434
00435 fDz = pDz ;
00436 fTthetaCphi = std::tan(pTheta)*std::cos(pPhi) ;
00437 fTthetaSphi = std::tan(pTheta)*std::sin(pPhi) ;
00438
00439 fDy1 = pDy ;
00440 fDx1 = pDx ;
00441 fDx2 = pDx ;
00442 fTalpha1 = std::tan(pAlpha) ;
00443
00444 fDy2 = pDy ;
00445 fDx3 = pDx ;
00446 fDx4 = pDx ;
00447 fTalpha2 = fTalpha1 ;
00448
00449 G4ThreeVector pt[8] ;
00450
00451 pt[0]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
00452 -fDz*fTthetaSphi-fDy1,-fDz);
00453 pt[1]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
00454 -fDz*fTthetaSphi-fDy1,-fDz);
00455 pt[2]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
00456 -fDz*fTthetaSphi+fDy1,-fDz);
00457 pt[3]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
00458 -fDz*fTthetaSphi+fDy1,-fDz);
00459 pt[4]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
00460 +fDz*fTthetaSphi-fDy2,+fDz);
00461 pt[5]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
00462 +fDz*fTthetaSphi-fDy2,+fDz);
00463 pt[6]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
00464 +fDz*fTthetaSphi+fDy2,+fDz);
00465 pt[7]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
00466 +fDz*fTthetaSphi+fDy2,+fDz);
00467
00468
00469
00470 good=MakePlane(pt[0],pt[4],pt[5],pt[1],fPlanes[0]);
00471 if (!good)
00472 {
00473 std::ostringstream message;
00474 message << "Face at ~-Y not planar for Solid: " << GetName();
00475 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
00476 FatalException, message);
00477 }
00478
00479
00480
00481 good=MakePlane(pt[2],pt[3],pt[7],pt[6],fPlanes[1]);
00482 if (!good)
00483 {
00484 std::ostringstream message;
00485 message << "Face at ~+Y not planar for Solid: " << GetName();
00486 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
00487 FatalException, message);
00488 }
00489
00490
00491
00492 good=MakePlane(pt[0],pt[2],pt[6],pt[4],fPlanes[2]);
00493 if (!good)
00494 {
00495 std::ostringstream message;
00496 message << "Face at ~-X not planar for Solid: " << GetName();
00497 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
00498 FatalException, message);
00499 }
00500
00501
00502
00503 good=MakePlane(pt[1],pt[5],pt[7],pt[3],fPlanes[3]);
00504 if (!good)
00505 {
00506 std::ostringstream message;
00507 message << "Face at ~+X not planar for Solid: " << GetName();
00508 G4Exception("G4Trap::G4Trap()", "GeomSolids0002",
00509 FatalException, message);
00510 }
00511 }
00512
00514
00515
00516
00517
00518
00519 G4Trap::G4Trap( const G4String& pName )
00520 : G4CSGSolid (pName), fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.),
00521 fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.),
00522 fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.)
00523 {
00524 MakePlanes();
00525 }
00526
00528
00529
00530
00531
00532 G4Trap::G4Trap( __void__& a )
00533 : G4CSGSolid(a), fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.),
00534 fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.),
00535 fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.)
00536 {
00537 MakePlanes();
00538 }
00539
00541
00542
00543
00544 G4Trap::~G4Trap()
00545 {
00546 }
00547
00549
00550
00551
00552 G4Trap::G4Trap(const G4Trap& rhs)
00553 : G4CSGSolid(rhs), fDz(rhs.fDz),
00554 fTthetaCphi(rhs.fTthetaCphi), fTthetaSphi(rhs.fTthetaSphi),
00555 fDy1(rhs.fDy1), fDx1(rhs.fDx1), fDx2(rhs.fDx2), fTalpha1(rhs.fTalpha1),
00556 fDy2(rhs.fDy2), fDx3(rhs.fDx3), fDx4(rhs.fDx4), fTalpha2(rhs.fTalpha2)
00557 {
00558 for (size_t i=0; i<4; ++i)
00559 {
00560 fPlanes[i].a = rhs.fPlanes[i].a;
00561 fPlanes[i].b = rhs.fPlanes[i].b;
00562 fPlanes[i].c = rhs.fPlanes[i].c;
00563 fPlanes[i].d = rhs.fPlanes[i].d;
00564 }
00565 }
00566
00568
00569
00570
00571 G4Trap& G4Trap::operator = (const G4Trap& rhs)
00572 {
00573
00574
00575 if (this == &rhs) { return *this; }
00576
00577
00578
00579 G4CSGSolid::operator=(rhs);
00580
00581
00582
00583 fDz = rhs.fDz;
00584 fTthetaCphi = rhs.fTthetaCphi; fTthetaSphi = rhs.fTthetaSphi;
00585 fDy1 = rhs.fDy1; fDx1 = rhs.fDx1; fDx2 = rhs.fDx2; fTalpha1 = rhs.fTalpha1;
00586 fDy2 = rhs.fDy2; fDx3 = rhs.fDx3; fDx4 = rhs.fDx4; fTalpha2 = rhs.fTalpha2;
00587 for (size_t i=0; i<4; ++i)
00588 {
00589 fPlanes[i].a = rhs.fPlanes[i].a;
00590 fPlanes[i].b = rhs.fPlanes[i].b;
00591 fPlanes[i].c = rhs.fPlanes[i].c;
00592 fPlanes[i].d = rhs.fPlanes[i].d;
00593 }
00594
00595 return *this;
00596 }
00597
00599
00600
00601
00602
00603 void G4Trap::SetAllParameters ( G4double pDz,
00604 G4double pTheta,
00605 G4double pPhi,
00606 G4double pDy1,
00607 G4double pDx1,
00608 G4double pDx2,
00609 G4double pAlp1,
00610 G4double pDy2,
00611 G4double pDx3,
00612 G4double pDx4,
00613 G4double pAlp2 )
00614 {
00615 if ( pDz<=0 || pDy1<=0 || pDx1<=0 || pDx2<=0 || pDy2<=0 || pDx3<=0 || pDx4<=0 )
00616 {
00617 std::ostringstream message;
00618 message << "Invalid Length Parameters for Solid: " << GetName() << G4endl
00619 << " X - "
00620 << pDx1 << ", " << pDx2 << ", " << pDx3 << ", " << pDx4 << G4endl
00621 << " Y - " << pDy1 << ", " << pDy2 << G4endl
00622 << " Z - " << pDz;
00623 G4Exception("G4Trap::SetAllParameters()", "GeomSolids0002",
00624 FatalException, message);
00625 }
00626 fCubicVolume= 0.;
00627 fSurfaceArea= 0.;
00628 fpPolyhedron = 0;
00629 fDz=pDz;
00630 fTthetaCphi=std::tan(pTheta)*std::cos(pPhi);
00631 fTthetaSphi=std::tan(pTheta)*std::sin(pPhi);
00632
00633 fDy1=pDy1;
00634 fDx1=pDx1;
00635 fDx2=pDx2;
00636 fTalpha1=std::tan(pAlp1);
00637
00638 fDy2=pDy2;
00639 fDx3=pDx3;
00640 fDx4=pDx4;
00641 fTalpha2=std::tan(pAlp2);
00642
00643 MakePlanes();
00644 }
00645
00647
00648
00649
00650 G4bool G4Trap::MakePlanes()
00651 {
00652 G4bool good = true;
00653
00654 G4ThreeVector pt[8] ;
00655
00656 pt[0]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
00657 -fDz*fTthetaSphi-fDy1,-fDz);
00658 pt[1]=G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
00659 -fDz*fTthetaSphi-fDy1,-fDz);
00660 pt[2]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
00661 -fDz*fTthetaSphi+fDy1,-fDz);
00662 pt[3]=G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
00663 -fDz*fTthetaSphi+fDy1,-fDz);
00664 pt[4]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
00665 +fDz*fTthetaSphi-fDy2,+fDz);
00666 pt[5]=G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
00667 +fDz*fTthetaSphi-fDy2,+fDz);
00668 pt[6]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
00669 +fDz*fTthetaSphi+fDy2,+fDz);
00670 pt[7]=G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
00671 +fDz*fTthetaSphi+fDy2,+fDz);
00672
00673
00674
00675 good=MakePlane(pt[0],pt[4],pt[5],pt[1],fPlanes[0]) ;
00676 if (!good)
00677 {
00678 std::ostringstream message;
00679 message << "Face at ~-Y not planar for Solid: " << GetName();
00680 G4Exception("G4Trap::MakePlanes()", "GeomSolids0002",
00681 FatalException, message);
00682 }
00683
00684
00685
00686 good=MakePlane(pt[2],pt[3],pt[7],pt[6],fPlanes[1]);
00687 if (!good)
00688 {
00689 std::ostringstream message;
00690 message << "Face at ~+Y not planar for Solid: " << GetName();
00691 G4Exception("G4Trap::MakePlanes()", "GeomSolids0002",
00692 FatalException, message);
00693 }
00694
00695
00696
00697 good=MakePlane(pt[0],pt[2],pt[6],pt[4],fPlanes[2]);
00698 if (!good)
00699 {
00700 std::ostringstream message;
00701 message << "Face at ~-X not planar for Solid: " << GetName();
00702 G4Exception("G4Trap::MakePlanes()", "GeomSolids0002",
00703 FatalException, message);
00704 }
00705
00706
00707
00708 good = MakePlane(pt[1],pt[5],pt[7],pt[3],fPlanes[3]);
00709 if ( !good )
00710 {
00711 std::ostringstream message;
00712 message << "Face at ~+X not planar for Solid: " << GetName();
00713 G4Exception("G4Trap::MakePlanes()", "GeomSolids0002",
00714 FatalException, message);
00715 }
00716
00717 return good;
00718 }
00719
00721
00722
00723
00724
00725
00726
00727
00728
00729 G4bool G4Trap::MakePlane( const G4ThreeVector& p1,
00730 const G4ThreeVector& p2,
00731 const G4ThreeVector& p3,
00732 const G4ThreeVector& p4,
00733 TrapSidePlane& plane )
00734 {
00735 G4double a, b, c, sd;
00736 G4ThreeVector v12, v13, v14, Vcross;
00737
00738 G4bool good;
00739
00740 v12 = p2 - p1;
00741 v13 = p3 - p1;
00742 v14 = p4 - p1;
00743 Vcross = v12.cross(v13);
00744
00745 if (std::fabs(Vcross.dot(v14)/(Vcross.mag()*v14.mag())) > kCoplanar_Tolerance)
00746 {
00747 good = false;
00748 }
00749 else
00750 {
00751
00752
00753
00754
00755
00756
00757
00758
00759
00760
00761
00762
00763
00764
00765 a = +(p4.y() - p2.y())*(p3.z() - p1.z())
00766 - (p3.y() - p1.y())*(p4.z() - p2.z());
00767
00768 b = -(p4.x() - p2.x())*(p3.z() - p1.z())
00769 + (p3.x() - p1.x())*(p4.z() - p2.z());
00770
00771 c = +(p4.x() - p2.x())*(p3.y() - p1.y())
00772 - (p3.x() - p1.x())*(p4.y() - p2.y());
00773
00774 sd = std::sqrt( a*a + b*b + c*c );
00775
00776 if( sd > 0 )
00777 {
00778 plane.a = a/sd;
00779 plane.b = b/sd;
00780 plane.c = c/sd;
00781 }
00782 else
00783 {
00784 std::ostringstream message;
00785 message << "Invalid parameters: norm.mod() <= 0, for Solid: "
00786 << GetName();
00787 G4Exception("G4Trap::MakePlanes()", "GeomSolids0002",
00788 FatalException, message) ;
00789 }
00790
00791
00792 plane.d = -( plane.a*p1.x() + plane.b*p1.y() + plane.c*p1.z() );
00793
00794 good = true;
00795 }
00796 return good;
00797 }
00798
00800
00801
00802
00803
00804 void G4Trap::ComputeDimensions( G4VPVParameterisation* p,
00805 const G4int n,
00806 const G4VPhysicalVolume* pRep )
00807 {
00808 p->ComputeDimensions(*this,n,pRep);
00809 }
00810
00811
00813
00814
00815
00816 G4bool G4Trap::CalculateExtent( const EAxis pAxis,
00817 const G4VoxelLimits& pVoxelLimit,
00818 const G4AffineTransform& pTransform,
00819 G4double& pMin, G4double& pMax) const
00820 {
00821 G4double xMin, xMax, yMin, yMax, zMin, zMax;
00822 G4bool flag;
00823
00824 if (!pTransform.IsRotated())
00825 {
00826
00827
00828
00829
00830 G4int i ;
00831 G4double xoffset;
00832 G4double yoffset;
00833 G4double zoffset;
00834 G4double temp[8] ;
00835 G4ThreeVector pt[8];
00836
00837 xoffset=pTransform.NetTranslation().x();
00838 yoffset=pTransform.NetTranslation().y();
00839 zoffset=pTransform.NetTranslation().z();
00840
00841 pt[0]=G4ThreeVector(xoffset-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
00842 yoffset-fDz*fTthetaSphi-fDy1,zoffset-fDz);
00843 pt[1]=G4ThreeVector(xoffset-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
00844 yoffset-fDz*fTthetaSphi-fDy1,zoffset-fDz);
00845 pt[2]=G4ThreeVector(xoffset-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
00846 yoffset-fDz*fTthetaSphi+fDy1,zoffset-fDz);
00847 pt[3]=G4ThreeVector(xoffset-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
00848 yoffset-fDz*fTthetaSphi+fDy1,zoffset-fDz);
00849 pt[4]=G4ThreeVector(xoffset+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
00850 yoffset+fDz*fTthetaSphi-fDy2,zoffset+fDz);
00851 pt[5]=G4ThreeVector(xoffset+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
00852 yoffset+fDz*fTthetaSphi-fDy2,zoffset+fDz);
00853 pt[6]=G4ThreeVector(xoffset+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
00854 yoffset+fDz*fTthetaSphi+fDy2,zoffset+fDz);
00855 pt[7]=G4ThreeVector(xoffset+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
00856 yoffset+fDz*fTthetaSphi+fDy2,zoffset+fDz);
00857 zMin=zoffset-fDz;
00858 zMax=zoffset+fDz;
00859
00860 if ( pVoxelLimit.IsZLimited() )
00861 {
00862 if ( (zMin > pVoxelLimit.GetMaxZExtent() + kCarTolerance)
00863 || (zMax < pVoxelLimit.GetMinZExtent() - kCarTolerance) )
00864 {
00865 return false;
00866 }
00867 else
00868 {
00869 if ( zMin < pVoxelLimit.GetMinZExtent() )
00870 {
00871 zMin = pVoxelLimit.GetMinZExtent() ;
00872 }
00873 if ( zMax > pVoxelLimit.GetMaxZExtent() )
00874 {
00875 zMax = pVoxelLimit.GetMaxZExtent() ;
00876 }
00877 }
00878 }
00879 temp[0] = pt[0].y()+(pt[4].y()-pt[0].y())*(zMin-pt[0].z())
00880 /(pt[4].z()-pt[0].z()) ;
00881 temp[1] = pt[0].y()+(pt[4].y()-pt[0].y())*(zMax-pt[0].z())
00882 /(pt[4].z()-pt[0].z()) ;
00883 temp[2] = pt[2].y()+(pt[6].y()-pt[2].y())*(zMin-pt[2].z())
00884 /(pt[6].z()-pt[2].z()) ;
00885 temp[3] = pt[2].y()+(pt[6].y()-pt[2].y())*(zMax-pt[2].z())
00886 /(pt[6].z()-pt[2].z()) ;
00887
00888 yMax = yoffset - std::fabs(fDz*fTthetaSphi) - fDy1 - fDy2 ;
00889 yMin = -yMax ;
00890
00891 for( i = 0 ; i < 4 ; i++ )
00892 {
00893 if( temp[i] > yMax ) yMax = temp[i] ;
00894 if( temp[i] < yMin ) yMin = temp[i] ;
00895 }
00896 if ( pVoxelLimit.IsYLimited() )
00897 {
00898 if ( (yMin > pVoxelLimit.GetMaxYExtent() + kCarTolerance)
00899 || (yMax < pVoxelLimit.GetMinYExtent() - kCarTolerance) )
00900 {
00901 return false;
00902 }
00903 else
00904 {
00905 if ( yMin < pVoxelLimit.GetMinYExtent() )
00906 {
00907 yMin = pVoxelLimit.GetMinYExtent() ;
00908 }
00909 if ( yMax > pVoxelLimit.GetMaxYExtent() )
00910 {
00911 yMax = pVoxelLimit.GetMaxYExtent() ;
00912 }
00913 }
00914 }
00915 temp[0] = pt[0].x()+(pt[4].x()-pt[0].x())
00916 *(zMin-pt[0].z())/(pt[4].z()-pt[0].z()) ;
00917 temp[1] = pt[0].x()+(pt[4].x()-pt[0].x())
00918 *(zMax-pt[0].z())/(pt[4].z()-pt[0].z()) ;
00919 temp[2] = pt[2].x()+(pt[6].x()-pt[2].x())
00920 *(zMin-pt[2].z())/(pt[6].z()-pt[2].z()) ;
00921 temp[3] = pt[2].x()+(pt[6].x()-pt[2].x())
00922 *(zMax-pt[2].z())/(pt[6].z()-pt[2].z()) ;
00923 temp[4] = pt[3].x()+(pt[7].x()-pt[3].x())
00924 *(zMin-pt[3].z())/(pt[7].z()-pt[3].z()) ;
00925 temp[5] = pt[3].x()+(pt[7].x()-pt[3].x())
00926 *(zMax-pt[3].z())/(pt[7].z()-pt[3].z()) ;
00927 temp[6] = pt[1].x()+(pt[5].x()-pt[1].x())
00928 *(zMin-pt[1].z())/(pt[5].z()-pt[1].z()) ;
00929 temp[7] = pt[1].x()+(pt[5].x()-pt[1].x())
00930 *(zMax-pt[1].z())/(pt[5].z()-pt[1].z()) ;
00931
00932 xMax = xoffset - std::fabs(fDz*fTthetaCphi) - fDx1 - fDx2 -fDx3 - fDx4 ;
00933 xMin = -xMax ;
00934
00935 for( i = 0 ; i < 8 ; i++ )
00936 {
00937 if( temp[i] > xMax) xMax = temp[i] ;
00938 if( temp[i] < xMin) xMin = temp[i] ;
00939 }
00940 if (pVoxelLimit.IsXLimited())
00941 {
00942 if ( (xMin > pVoxelLimit.GetMaxXExtent() + kCarTolerance)
00943 || (xMax < pVoxelLimit.GetMinXExtent() - kCarTolerance) )
00944 {
00945 return false;
00946 }
00947 else
00948 {
00949 if ( xMin < pVoxelLimit.GetMinXExtent() )
00950 {
00951 xMin = pVoxelLimit.GetMinXExtent() ;
00952 }
00953 if ( xMax > pVoxelLimit.GetMaxXExtent() )
00954 {
00955 xMax = pVoxelLimit.GetMaxXExtent() ;
00956 }
00957 }
00958 }
00959 switch (pAxis)
00960 {
00961 case kXAxis:
00962 pMin=xMin;
00963 pMax=xMax;
00964 break;
00965
00966 case kYAxis:
00967 pMin=yMin;
00968 pMax=yMax;
00969 break;
00970
00971 case kZAxis:
00972 pMin=zMin;
00973 pMax=zMax;
00974 break;
00975
00976 default:
00977 break;
00978 }
00979 pMin -= kCarTolerance;
00980 pMax += kCarTolerance;
00981
00982 flag = true;
00983 }
00984 else
00985 {
00986 G4bool existsAfterClip = false ;
00987 G4ThreeVectorList* vertices;
00988 pMin = +kInfinity;
00989 pMax = -kInfinity;
00990
00991
00992
00993 vertices = CreateRotatedVertices(pTransform);
00994
00995 xMin = +kInfinity; yMin = +kInfinity; zMin = +kInfinity;
00996 xMax = -kInfinity; yMax = -kInfinity; zMax = -kInfinity;
00997
00998 for( G4int nv = 0 ; nv < 8 ; nv++ )
00999 {
01000 if( (*vertices)[nv].x() > xMax ) xMax = (*vertices)[nv].x();
01001 if( (*vertices)[nv].y() > yMax ) yMax = (*vertices)[nv].y();
01002 if( (*vertices)[nv].z() > zMax ) zMax = (*vertices)[nv].z();
01003
01004 if( (*vertices)[nv].x() < xMin ) xMin = (*vertices)[nv].x();
01005 if( (*vertices)[nv].y() < yMin ) yMin = (*vertices)[nv].y();
01006 if( (*vertices)[nv].z() < zMin ) zMin = (*vertices)[nv].z();
01007 }
01008 if ( pVoxelLimit.IsZLimited() )
01009 {
01010 if ( (zMin > pVoxelLimit.GetMaxZExtent() + kCarTolerance)
01011 || (zMax < pVoxelLimit.GetMinZExtent() - kCarTolerance) )
01012 {
01013 delete vertices ;
01014 return false;
01015 }
01016 else
01017 {
01018 if ( zMin < pVoxelLimit.GetMinZExtent() )
01019 {
01020 zMin = pVoxelLimit.GetMinZExtent() ;
01021 }
01022 if ( zMax > pVoxelLimit.GetMaxZExtent() )
01023 {
01024 zMax = pVoxelLimit.GetMaxZExtent() ;
01025 }
01026 }
01027 }
01028 if ( pVoxelLimit.IsYLimited() )
01029 {
01030 if ( (yMin > pVoxelLimit.GetMaxYExtent() + kCarTolerance)
01031 || (yMax < pVoxelLimit.GetMinYExtent() - kCarTolerance) )
01032 {
01033 delete vertices ;
01034 return false;
01035 }
01036 else
01037 {
01038 if ( yMin < pVoxelLimit.GetMinYExtent() )
01039 {
01040 yMin = pVoxelLimit.GetMinYExtent() ;
01041 }
01042 if ( yMax > pVoxelLimit.GetMaxYExtent() )
01043 {
01044 yMax = pVoxelLimit.GetMaxYExtent() ;
01045 }
01046 }
01047 }
01048 if ( pVoxelLimit.IsXLimited() )
01049 {
01050 if ( (xMin > pVoxelLimit.GetMaxXExtent() + kCarTolerance)
01051 || (xMax < pVoxelLimit.GetMinXExtent() - kCarTolerance) )
01052 {
01053 delete vertices ;
01054 return false ;
01055 }
01056 else
01057 {
01058 if ( xMin < pVoxelLimit.GetMinXExtent() )
01059 {
01060 xMin = pVoxelLimit.GetMinXExtent() ;
01061 }
01062 if ( xMax > pVoxelLimit.GetMaxXExtent() )
01063 {
01064 xMax = pVoxelLimit.GetMaxXExtent() ;
01065 }
01066 }
01067 }
01068 switch (pAxis)
01069 {
01070 case kXAxis:
01071 pMin=xMin;
01072 pMax=xMax;
01073 break;
01074
01075 case kYAxis:
01076 pMin=yMin;
01077 pMax=yMax;
01078 break;
01079
01080 case kZAxis:
01081 pMin=zMin;
01082 pMax=zMax;
01083 break;
01084
01085 default:
01086 break;
01087 }
01088 if ( (pMin != kInfinity) || (pMax != -kInfinity) )
01089 {
01090 existsAfterClip=true;
01091
01092
01093
01094 pMin -= kCarTolerance ;
01095 pMax += kCarTolerance ;
01096 }
01097 delete vertices ;
01098 flag = existsAfterClip ;
01099 }
01100 return flag;
01101 }
01102
01103
01105
01106
01107
01108 EInside G4Trap::Inside( const G4ThreeVector& p ) const
01109 {
01110 EInside in;
01111 G4double Dist;
01112 G4int i;
01113 if ( std::fabs(p.z()) <= fDz-kCarTolerance*0.5)
01114 {
01115 in = kInside;
01116
01117 for ( i = 0;i < 4;i++ )
01118 {
01119 Dist = fPlanes[i].a*p.x() + fPlanes[i].b*p.y()
01120 +fPlanes[i].c*p.z() + fPlanes[i].d;
01121
01122 if (Dist > kCarTolerance*0.5) return in = kOutside;
01123 else if (Dist > -kCarTolerance*0.5) in = kSurface;
01124
01125 }
01126 }
01127 else if (std::fabs(p.z()) <= fDz+kCarTolerance*0.5)
01128 {
01129 in = kSurface;
01130
01131 for ( i = 0; i < 4; i++ )
01132 {
01133 Dist = fPlanes[i].a*p.x() + fPlanes[i].b*p.y()
01134 +fPlanes[i].c*p.z() + fPlanes[i].d;
01135
01136 if (Dist > kCarTolerance*0.5) return in = kOutside;
01137 }
01138 }
01139 else in = kOutside;
01140
01141 return in;
01142 }
01143
01145
01146
01147
01148
01149 G4ThreeVector G4Trap::SurfaceNormal( const G4ThreeVector& p ) const
01150 {
01151 G4int i, noSurfaces = 0;
01152 G4double dist, distz, distx, disty, distmx, distmy, safe = kInfinity;
01153 G4double delta = 0.5*kCarTolerance;
01154 G4ThreeVector norm, sumnorm(0.,0.,0.);
01155
01156 for (i = 0; i < 4; i++)
01157 {
01158 dist = std::fabs(fPlanes[i].a*p.x() + fPlanes[i].b*p.y()
01159 + fPlanes[i].c*p.z() + fPlanes[i].d);
01160 if ( dist < safe )
01161 {
01162 safe = dist;
01163 }
01164 }
01165 distz = std::fabs( std::fabs( p.z() ) - fDz );
01166
01167 distmy = std::fabs( fPlanes[0].a*p.x() + fPlanes[0].b*p.y()
01168 + fPlanes[0].c*p.z() + fPlanes[0].d );
01169
01170 disty = std::fabs( fPlanes[1].a*p.x() + fPlanes[1].b*p.y()
01171 + fPlanes[1].c*p.z() + fPlanes[1].d );
01172
01173 distmx = std::fabs( fPlanes[2].a*p.x() + fPlanes[2].b*p.y()
01174 + fPlanes[2].c*p.z() + fPlanes[2].d );
01175
01176 distx = std::fabs( fPlanes[3].a*p.x() + fPlanes[3].b*p.y()
01177 + fPlanes[3].c*p.z() + fPlanes[3].d );
01178
01179 G4ThreeVector nX = G4ThreeVector(fPlanes[3].a,fPlanes[3].b,fPlanes[3].c);
01180 G4ThreeVector nmX = G4ThreeVector(fPlanes[2].a,fPlanes[2].b,fPlanes[2].c);
01181 G4ThreeVector nY = G4ThreeVector(fPlanes[1].a,fPlanes[1].b,fPlanes[1].c);
01182 G4ThreeVector nmY = G4ThreeVector(fPlanes[0].a,fPlanes[0].b,fPlanes[0].c);
01183 G4ThreeVector nZ = G4ThreeVector(0.,0.,1.0);
01184
01185 if (distx <= delta)
01186 {
01187 noSurfaces ++;
01188 sumnorm += nX;
01189 }
01190 if (distmx <= delta)
01191 {
01192 noSurfaces ++;
01193 sumnorm += nmX;
01194 }
01195 if (disty <= delta)
01196 {
01197 noSurfaces ++;
01198 sumnorm += nY;
01199 }
01200 if (distmy <= delta)
01201 {
01202 noSurfaces ++;
01203 sumnorm += nmY;
01204 }
01205 if (distz <= delta)
01206 {
01207 noSurfaces ++;
01208 if ( p.z() >= 0.) sumnorm += nZ;
01209 else sumnorm -= nZ;
01210 }
01211 if ( noSurfaces == 0 )
01212 {
01213 #ifdef G4CSGDEBUG
01214 G4Exception("G4Trap::SurfaceNormal(p)", "GeomSolids1002",
01215 JustWarning, "Point p is not on surface !?" );
01216 #endif
01217 norm = ApproxSurfaceNormal(p);
01218 }
01219 else if ( noSurfaces == 1 ) norm = sumnorm;
01220 else norm = sumnorm.unit();
01221 return norm;
01222 }
01223
01225
01226
01227
01228
01229 G4ThreeVector G4Trap::ApproxSurfaceNormal( const G4ThreeVector& p ) const
01230 {
01231 G4double safe=kInfinity,Dist,safez;
01232 G4int i,imin=0;
01233 for (i=0;i<4;i++)
01234 {
01235 Dist=std::fabs(fPlanes[i].a*p.x()+fPlanes[i].b*p.y()
01236 +fPlanes[i].c*p.z()+fPlanes[i].d);
01237 if (Dist<safe)
01238 {
01239 safe=Dist;
01240 imin=i;
01241 }
01242 }
01243 safez=std::fabs(std::fabs(p.z())-fDz);
01244 if (safe<safez)
01245 {
01246 return G4ThreeVector(fPlanes[imin].a,fPlanes[imin].b,fPlanes[imin].c);
01247 }
01248 else
01249 {
01250 if (p.z()>0)
01251 {
01252 return G4ThreeVector(0,0,1);
01253 }
01254 else
01255 {
01256 return G4ThreeVector(0,0,-1);
01257 }
01258 }
01259 }
01260
01262
01263
01264
01265
01266
01267
01268
01269
01270 G4double G4Trap::DistanceToIn( const G4ThreeVector& p,
01271 const G4ThreeVector& v ) const
01272 {
01273
01274 G4double snxt;
01275 G4double max,smax,smin;
01276 G4double pdist,Comp,vdist;
01277 G4int i;
01278
01279
01280
01281 if ( v.z() > 0 )
01282 {
01283 max = fDz - p.z() ;
01284 if (max > 0.5*kCarTolerance)
01285 {
01286 smax = max/v.z();
01287 smin = (-fDz-p.z())/v.z();
01288 }
01289 else
01290 {
01291 return snxt=kInfinity;
01292 }
01293 }
01294 else if (v.z() < 0 )
01295 {
01296 max = - fDz - p.z() ;
01297 if (max < -0.5*kCarTolerance )
01298 {
01299 smax=max/v.z();
01300 smin=(fDz-p.z())/v.z();
01301 }
01302 else
01303 {
01304 return snxt=kInfinity;
01305 }
01306 }
01307 else
01308 {
01309 if (std::fabs(p.z())<fDz - 0.5*kCarTolerance)
01310 {
01311 smin=0;
01312 smax=kInfinity;
01313 }
01314 else
01315 {
01316 return snxt=kInfinity;
01317 }
01318 }
01319
01320 for (i=0;i<4;i++)
01321 {
01322 pdist=fPlanes[i].a*p.x()+fPlanes[i].b*p.y()
01323 +fPlanes[i].c*p.z()+fPlanes[i].d;
01324 Comp=fPlanes[i].a*v.x()+fPlanes[i].b*v.y()+fPlanes[i].c*v.z();
01325 if ( pdist >= -0.5*kCarTolerance )
01326 {
01327
01328
01329
01330 if (Comp >= 0)
01331 {
01332 return snxt=kInfinity ;
01333 }
01334 else
01335 {
01336 vdist=-pdist/Comp;
01337 if (vdist>smin)
01338 {
01339 if (vdist<smax)
01340 {
01341 smin = vdist;
01342 }
01343 else
01344 {
01345 return snxt=kInfinity;
01346 }
01347 }
01348 }
01349 }
01350 else
01351 {
01352
01353
01354
01355 if (Comp>0)
01356 {
01357 vdist=-pdist/Comp;
01358 if (vdist<smax)
01359 {
01360 if (vdist>smin)
01361 {
01362 smax=vdist;
01363 }
01364 else
01365 {
01366 return snxt=kInfinity;
01367 }
01368 }
01369 }
01370 }
01371 }
01372
01373
01374
01375 if (smin >=0 )
01376 {
01377 snxt = smin ;
01378 }
01379 else
01380 {
01381 snxt = 0 ;
01382 }
01383 return snxt;
01384 }
01385
01387
01388
01389
01390
01391
01392 G4double G4Trap::DistanceToIn( const G4ThreeVector& p ) const
01393 {
01394 G4double safe=0.0,Dist;
01395 G4int i;
01396 safe=std::fabs(p.z())-fDz;
01397 for (i=0;i<4;i++)
01398 {
01399 Dist=fPlanes[i].a*p.x()+fPlanes[i].b*p.y()
01400 +fPlanes[i].c*p.z()+fPlanes[i].d;
01401 if (Dist > safe) safe=Dist;
01402 }
01403 if (safe<0) safe=0;
01404 return safe;
01405 }
01406
01408
01409
01410
01411
01412 G4double G4Trap::DistanceToOut(const G4ThreeVector& p, const G4ThreeVector& v,
01413 const G4bool calcNorm,
01414 G4bool *validNorm, G4ThreeVector *n) const
01415 {
01416 Eside side = kUndef;
01417 G4double snxt;
01418 G4double pdist,Comp,vdist,max;
01419
01420
01421
01422 if (v.z()>0)
01423 {
01424 max=fDz-p.z();
01425 if (max>kCarTolerance/2)
01426 {
01427 snxt=max/v.z();
01428 side=kPZ;
01429 }
01430 else
01431 {
01432 if (calcNorm)
01433 {
01434 *validNorm=true;
01435 *n=G4ThreeVector(0,0,1);
01436 }
01437 return snxt=0;
01438 }
01439 }
01440 else if (v.z()<0)
01441 {
01442 max=-fDz-p.z();
01443 if (max<-kCarTolerance/2)
01444 {
01445 snxt=max/v.z();
01446 side=kMZ;
01447 }
01448 else
01449 {
01450 if (calcNorm)
01451 {
01452 *validNorm=true;
01453 *n=G4ThreeVector(0,0,-1);
01454 }
01455 return snxt=0;
01456 }
01457 }
01458 else
01459 {
01460 snxt=kInfinity;
01461 }
01462
01463
01464
01465
01466 pdist=fPlanes[0].a*p.x()+fPlanes[0].b*p.y()+fPlanes[0].c*p.z()+fPlanes[0].d;
01467 Comp=fPlanes[0].a*v.x()+fPlanes[0].b*v.y()+fPlanes[0].c*v.z();
01468 if (pdist>0)
01469 {
01470
01471 if (Comp>0)
01472 {
01473
01474 if (calcNorm)
01475 {
01476 *validNorm=true;
01477 *n=G4ThreeVector(fPlanes[0].a,fPlanes[0].b,fPlanes[0].c);
01478 }
01479 return snxt=0;
01480 }
01481 }
01482 else if (pdist<-kCarTolerance/2)
01483 {
01484
01485 if (Comp>0)
01486 {
01487
01488 vdist=-pdist/Comp;
01489 if (vdist<snxt)
01490 {
01491 snxt=vdist;
01492 side=ks0;
01493 }
01494 }
01495 }
01496 else
01497 {
01498
01499 if (Comp>0)
01500 {
01501 if (calcNorm)
01502 {
01503 *validNorm=true;
01504 *n=G4ThreeVector(fPlanes[0].a,fPlanes[0].b,fPlanes[0].c);
01505 }
01506 return snxt=0;
01507 }
01508 }
01509
01510
01511
01512
01513 pdist=fPlanes[1].a*p.x()+fPlanes[1].b*p.y()+fPlanes[1].c*p.z()+fPlanes[1].d;
01514 Comp=fPlanes[1].a*v.x()+fPlanes[1].b*v.y()+fPlanes[1].c*v.z();
01515 if (pdist>0)
01516 {
01517
01518 if (Comp>0)
01519 {
01520
01521 if (calcNorm)
01522 {
01523 *validNorm=true;
01524 *n=G4ThreeVector(fPlanes[1].a,fPlanes[1].b,fPlanes[1].c);
01525 }
01526 return snxt=0;
01527 }
01528 }
01529 else if (pdist<-kCarTolerance/2)
01530 {
01531
01532 if (Comp>0)
01533 {
01534
01535 vdist=-pdist/Comp;
01536 if (vdist<snxt)
01537 {
01538 snxt=vdist;
01539 side=ks1;
01540 }
01541 }
01542 }
01543 else
01544 {
01545
01546 if (Comp>0)
01547 {
01548 if (calcNorm)
01549 {
01550 *validNorm=true;
01551 *n=G4ThreeVector(fPlanes[1].a,fPlanes[1].b,fPlanes[1].c);
01552 }
01553 return snxt=0;
01554 }
01555 }
01556
01557
01558
01559
01560 pdist=fPlanes[2].a*p.x()+fPlanes[2].b*p.y()+fPlanes[2].c*p.z()+fPlanes[2].d;
01561 Comp=fPlanes[2].a*v.x()+fPlanes[2].b*v.y()+fPlanes[2].c*v.z();
01562 if (pdist>0)
01563 {
01564
01565 if (Comp>0)
01566 {
01567
01568 if (calcNorm)
01569 {
01570 *validNorm=true;
01571 *n=G4ThreeVector(fPlanes[2].a,fPlanes[2].b,fPlanes[2].c);
01572 }
01573 return snxt=0;
01574 }
01575 }
01576 else if (pdist<-kCarTolerance/2)
01577 {
01578
01579 if (Comp>0)
01580 {
01581
01582 vdist=-pdist/Comp;
01583 if (vdist<snxt)
01584 {
01585 snxt=vdist;
01586 side=ks2;
01587 }
01588 }
01589 }
01590 else
01591 {
01592
01593 if (Comp>0)
01594 {
01595 if (calcNorm)
01596 {
01597 *validNorm=true;
01598 *n=G4ThreeVector(fPlanes[2].a,fPlanes[2].b,fPlanes[2].c);
01599 }
01600 return snxt=0;
01601 }
01602 }
01603
01604
01605
01606
01607 pdist=fPlanes[3].a*p.x()+fPlanes[3].b*p.y()+fPlanes[3].c*p.z()+fPlanes[3].d;
01608 Comp=fPlanes[3].a*v.x()+fPlanes[3].b*v.y()+fPlanes[3].c*v.z();
01609 if (pdist>0)
01610 {
01611
01612 if (Comp>0)
01613 {
01614
01615 if (calcNorm)
01616 {
01617 *validNorm=true;
01618 *n=G4ThreeVector(fPlanes[3].a,fPlanes[3].b,fPlanes[3].c);
01619 }
01620 return snxt=0;
01621 }
01622 }
01623 else if (pdist<-kCarTolerance/2)
01624 {
01625
01626 if (Comp>0)
01627 {
01628
01629 vdist=-pdist/Comp;
01630 if (vdist<snxt)
01631 {
01632 snxt=vdist;
01633 side=ks3;
01634 }
01635 }
01636 }
01637 else
01638 {
01639
01640 if (Comp>0)
01641 {
01642 if (calcNorm)
01643 {
01644 *validNorm=true;
01645 *n=G4ThreeVector(fPlanes[3].a,fPlanes[3].b,fPlanes[3].c);
01646 }
01647 return snxt=0;
01648 }
01649 }
01650
01651
01652 if (calcNorm)
01653 {
01654 *validNorm=true;
01655 switch(side)
01656 {
01657 case ks0:
01658 *n=G4ThreeVector(fPlanes[0].a,fPlanes[0].b,fPlanes[0].c);
01659 break;
01660 case ks1:
01661 *n=G4ThreeVector(fPlanes[1].a,fPlanes[1].b,fPlanes[1].c);
01662 break;
01663 case ks2:
01664 *n=G4ThreeVector(fPlanes[2].a,fPlanes[2].b,fPlanes[2].c);
01665 break;
01666 case ks3:
01667 *n=G4ThreeVector(fPlanes[3].a,fPlanes[3].b,fPlanes[3].c);
01668 break;
01669 case kMZ:
01670 *n=G4ThreeVector(0,0,-1);
01671 break;
01672 case kPZ:
01673 *n=G4ThreeVector(0,0,1);
01674 break;
01675 default:
01676 G4cout << G4endl;
01677 DumpInfo();
01678 std::ostringstream message;
01679 G4int oldprc = message.precision(16);
01680 message << "Undefined side for valid surface normal to solid."
01681 << G4endl
01682 << "Position:" << G4endl << G4endl
01683 << "p.x() = " << p.x()/mm << " mm" << G4endl
01684 << "p.y() = " << p.y()/mm << " mm" << G4endl
01685 << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
01686 << "Direction:" << G4endl << G4endl
01687 << "v.x() = " << v.x() << G4endl
01688 << "v.y() = " << v.y() << G4endl
01689 << "v.z() = " << v.z() << G4endl << G4endl
01690 << "Proposed distance :" << G4endl << G4endl
01691 << "snxt = " << snxt/mm << " mm" << G4endl;
01692 message.precision(oldprc);
01693 G4Exception("G4Trap::DistanceToOut(p,v,..)","GeomSolids1002",
01694 JustWarning, message);
01695 break;
01696 }
01697 }
01698 return snxt;
01699 }
01700
01702
01703
01704
01705
01706 G4double G4Trap::DistanceToOut( const G4ThreeVector& p ) const
01707 {
01708 G4double safe=0.0,Dist;
01709 G4int i;
01710
01711 #ifdef G4CSGDEBUG
01712 if( Inside(p) == kOutside )
01713 {
01714 G4int oldprc = G4cout.precision(16) ;
01715 G4cout << G4endl ;
01716 DumpInfo();
01717 G4cout << "Position:" << G4endl << G4endl ;
01718 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
01719 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
01720 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
01721 G4cout.precision(oldprc) ;
01722 G4Exception("G4Trap::DistanceToOut(p)",
01723 "GeomSolids1002", JustWarning, "Point p is outside !?" );
01724 }
01725 #endif
01726
01727 safe=fDz-std::fabs(p.z());
01728 if (safe<0) safe=0;
01729 else
01730 {
01731 for (i=0;i<4;i++)
01732 {
01733 Dist=-(fPlanes[i].a*p.x()+fPlanes[i].b*p.y()
01734 +fPlanes[i].c*p.z()+fPlanes[i].d);
01735 if (Dist<safe) safe=Dist;
01736 }
01737 if (safe<0) safe=0;
01738 }
01739 return safe;
01740 }
01741
01743
01744
01745
01746
01747
01748
01749
01750
01751 G4ThreeVectorList*
01752 G4Trap::CreateRotatedVertices( const G4AffineTransform& pTransform ) const
01753 {
01754 G4ThreeVectorList *vertices;
01755 vertices=new G4ThreeVectorList();
01756 if (vertices)
01757 {
01758 vertices->reserve(8);
01759 G4ThreeVector vertex0(-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
01760 -fDz*fTthetaSphi-fDy1,-fDz);
01761 G4ThreeVector vertex1(-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
01762 -fDz*fTthetaSphi-fDy1,-fDz);
01763 G4ThreeVector vertex2(-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
01764 -fDz*fTthetaSphi+fDy1,-fDz);
01765 G4ThreeVector vertex3(-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
01766 -fDz*fTthetaSphi+fDy1,-fDz);
01767 G4ThreeVector vertex4(+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
01768 +fDz*fTthetaSphi-fDy2,+fDz);
01769 G4ThreeVector vertex5(+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
01770 +fDz*fTthetaSphi-fDy2,+fDz);
01771 G4ThreeVector vertex6(+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
01772 +fDz*fTthetaSphi+fDy2,+fDz);
01773 G4ThreeVector vertex7(+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
01774 +fDz*fTthetaSphi+fDy2,+fDz);
01775
01776 vertices->push_back(pTransform.TransformPoint(vertex0));
01777 vertices->push_back(pTransform.TransformPoint(vertex1));
01778 vertices->push_back(pTransform.TransformPoint(vertex2));
01779 vertices->push_back(pTransform.TransformPoint(vertex3));
01780 vertices->push_back(pTransform.TransformPoint(vertex4));
01781 vertices->push_back(pTransform.TransformPoint(vertex5));
01782 vertices->push_back(pTransform.TransformPoint(vertex6));
01783 vertices->push_back(pTransform.TransformPoint(vertex7));
01784 }
01785 else
01786 {
01787 DumpInfo();
01788 G4Exception("G4Trap::CreateRotatedVertices()",
01789 "GeomSolids0003", FatalException,
01790 "Error in allocation of vertices. Out of memory !");
01791 }
01792 return vertices;
01793 }
01794
01796
01797
01798
01799 G4GeometryType G4Trap::GetEntityType() const
01800 {
01801 return G4String("G4Trap");
01802 }
01803
01805
01806
01807
01808 G4VSolid* G4Trap::Clone() const
01809 {
01810 return new G4Trap(*this);
01811 }
01812
01814
01815
01816
01817 std::ostream& G4Trap::StreamInfo( std::ostream& os ) const
01818 {
01819 G4int oldprc = os.precision(16);
01820 os << "-----------------------------------------------------------\n"
01821 << " *** Dump for solid - " << GetName() << " ***\n"
01822 << " ===================================================\n"
01823 << " Solid type: G4Trap\n"
01824 << " Parameters: \n"
01825 << " half length Z: " << fDz/mm << " mm \n"
01826 << " half length Y of face -fDz: " << fDy1/mm << " mm \n"
01827 << " half length X of side -fDy1, face -fDz: " << fDx1/mm << " mm \n"
01828 << " half length X of side +fDy1, face -fDz: " << fDx2/mm << " mm \n"
01829 << " half length Y of face +fDz: " << fDy2/mm << " mm \n"
01830 << " half length X of side -fDy2, face +fDz: " << fDx3/mm << " mm \n"
01831 << " half length X of side +fDy2, face +fDz: " << fDx4/mm << " mm \n"
01832 << " std::tan(theta)*std::cos(phi): " << fTthetaCphi/degree << " degrees \n"
01833 << " std::tan(theta)*std::sin(phi): " << fTthetaSphi/degree << " degrees \n"
01834 << " std::tan(alpha), -fDz: " << fTalpha1/degree << " degrees \n"
01835 << " std::tan(alpha), +fDz: " << fTalpha2/degree << " degrees \n"
01836 << " trap side plane equations:\n"
01837 << " " << fPlanes[0].a << " X + " << fPlanes[0].b << " Y + "
01838 << fPlanes[0].c << " Z + " << fPlanes[0].d << " = 0\n"
01839 << " " << fPlanes[1].a << " X + " << fPlanes[1].b << " Y + "
01840 << fPlanes[1].c << " Z + " << fPlanes[1].d << " = 0\n"
01841 << " " << fPlanes[2].a << " X + " << fPlanes[2].b << " Y + "
01842 << fPlanes[2].c << " Z + " << fPlanes[2].d << " = 0\n"
01843 << " " << fPlanes[3].a << " X + " << fPlanes[3].b << " Y + "
01844 << fPlanes[3].c << " Z + " << fPlanes[3].d << " = 0\n"
01845 << "-----------------------------------------------------------\n";
01846 os.precision(oldprc);
01847
01848 return os;
01849 }
01850
01852
01853
01854
01855
01856
01857 G4ThreeVector G4Trap::GetPointOnPlane(G4ThreeVector p0, G4ThreeVector p1,
01858 G4ThreeVector p2, G4ThreeVector p3,
01859 G4double& area) const
01860 {
01861 G4double lambda1, lambda2, chose, aOne, aTwo;
01862 G4ThreeVector t, u, v, w, Area, normal;
01863
01864 t = p1 - p0;
01865 u = p2 - p1;
01866 v = p3 - p2;
01867 w = p0 - p3;
01868
01869 Area = G4ThreeVector(w.y()*v.z() - w.z()*v.y(),
01870 w.z()*v.x() - w.x()*v.z(),
01871 w.x()*v.y() - w.y()*v.x());
01872
01873 aOne = 0.5*Area.mag();
01874
01875 Area = G4ThreeVector(t.y()*u.z() - t.z()*u.y(),
01876 t.z()*u.x() - t.x()*u.z(),
01877 t.x()*u.y() - t.y()*u.x());
01878
01879 aTwo = 0.5*Area.mag();
01880
01881 area = aOne + aTwo;
01882
01883 chose = RandFlat::shoot(0.,aOne+aTwo);
01884
01885 if( (chose>=0.) && (chose < aOne) )
01886 {
01887 lambda1 = RandFlat::shoot(0.,1.);
01888 lambda2 = RandFlat::shoot(0.,lambda1);
01889 return (p2+lambda1*v+lambda2*w);
01890 }
01891
01892
01893
01894 lambda1 = RandFlat::shoot(0.,1.);
01895 lambda2 = RandFlat::shoot(0.,lambda1);
01896
01897 return (p0+lambda1*t+lambda2*u);
01898 }
01899
01901
01902
01903
01904 G4ThreeVector G4Trap::GetPointOnSurface() const
01905 {
01906 G4double aOne, aTwo, aThree, aFour, aFive, aSix, chose;
01907 G4ThreeVector One, Two, Three, Four, Five, Six, test;
01908 G4ThreeVector pt[8];
01909
01910 pt[0] = G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1-fDx1,
01911 -fDz*fTthetaSphi-fDy1,-fDz);
01912 pt[1] = G4ThreeVector(-fDz*fTthetaCphi-fDy1*fTalpha1+fDx1,
01913 -fDz*fTthetaSphi-fDy1,-fDz);
01914 pt[2] = G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1-fDx2,
01915 -fDz*fTthetaSphi+fDy1,-fDz);
01916 pt[3] = G4ThreeVector(-fDz*fTthetaCphi+fDy1*fTalpha1+fDx2,
01917 -fDz*fTthetaSphi+fDy1,-fDz);
01918 pt[4] = G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2-fDx3,
01919 +fDz*fTthetaSphi-fDy2,+fDz);
01920 pt[5] = G4ThreeVector(+fDz*fTthetaCphi-fDy2*fTalpha2+fDx3,
01921 +fDz*fTthetaSphi-fDy2,+fDz);
01922 pt[6] = G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2-fDx4,
01923 +fDz*fTthetaSphi+fDy2,+fDz);
01924 pt[7] = G4ThreeVector(+fDz*fTthetaCphi+fDy2*fTalpha2+fDx4,
01925 +fDz*fTthetaSphi+fDy2,+fDz);
01926
01927
01928
01929 One = GetPointOnPlane(pt[0],pt[1],pt[3],pt[2], aOne);
01930 Two = GetPointOnPlane(pt[4],pt[5],pt[7],pt[6], aTwo);
01931 Three = GetPointOnPlane(pt[6],pt[7],pt[3],pt[2], aThree);
01932 Four = GetPointOnPlane(pt[4],pt[5],pt[1],pt[0], aFour);
01933 Five = GetPointOnPlane(pt[0],pt[2],pt[6],pt[4], aFive);
01934 Six = GetPointOnPlane(pt[1],pt[3],pt[7],pt[5], aSix);
01935
01936 chose = RandFlat::shoot(0.,aOne+aTwo+aThree+aFour+aFive+aSix);
01937 if( (chose>=0.) && (chose<aOne) )
01938 { return One; }
01939 else if( (chose>=aOne) && (chose<aOne+aTwo) )
01940 { return Two; }
01941 else if( (chose>=aOne+aTwo) && (chose<aOne+aTwo+aThree) )
01942 { return Three; }
01943 else if( (chose>=aOne+aTwo+aThree) && (chose<aOne+aTwo+aThree+aFour) )
01944 { return Four; }
01945 else if( (chose>=aOne+aTwo+aThree+aFour)
01946 && (chose<aOne+aTwo+aThree+aFour+aFive) )
01947 { return Five; }
01948 return Six;
01949 }
01950
01952
01953
01954
01955 void G4Trap::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
01956 {
01957 scene.AddSolid (*this);
01958 }
01959
01960 G4Polyhedron* G4Trap::CreatePolyhedron () const
01961 {
01962 G4double phi = std::atan2(fTthetaSphi, fTthetaCphi);
01963 G4double alpha1 = std::atan(fTalpha1);
01964 G4double alpha2 = std::atan(fTalpha2);
01965 G4double theta = std::atan(std::sqrt(fTthetaCphi*fTthetaCphi+fTthetaSphi*fTthetaSphi));
01966
01967 return new G4PolyhedronTrap(fDz, theta, phi,
01968 fDy1, fDx1, fDx2, alpha1,
01969 fDy2, fDx3, fDx4, alpha2);
01970 }
01971
01972 G4NURBS* G4Trap::CreateNURBS () const
01973 {
01974
01975 return 0 ;
01976 }