G4TwistTubsFlatSide.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // $Id: G4TwistTubsFlatSide.cc 67011 2013-01-29 16:17:41Z gcosmo $
00028 //
00029 // 
00030 // --------------------------------------------------------------------
00031 // GEANT 4 class source file
00032 //
00033 //
00034 // G4TwistTubsFlatSide.cc
00035 //
00036 // Author: 
00037 //   01-Aug-2002 - Kotoyo Hoshina (hoshina@hepburn.s.chiba-u.ac.jp)
00038 //
00039 // History:
00040 //   13-Nov-2003 - O.Link (Oliver.Link@cern.ch), Integration in Geant4
00041 //                 from original version in Jupiter-2.5.02 application.
00042 // --------------------------------------------------------------------
00043 
00044 #include "G4TwistTubsFlatSide.hh"
00045 #include "G4GeometryTolerance.hh"
00046 
00047 //=====================================================================
00048 //* constructors ------------------------------------------------------
00049 
00050 G4TwistTubsFlatSide::G4TwistTubsFlatSide(const G4String         &name,
00051                              const G4RotationMatrix &rot,
00052                              const G4ThreeVector    &tlate,
00053                              const G4ThreeVector    &n,
00054                              const EAxis             axis0 ,
00055                              const EAxis             axis1 ,
00056                                    G4double          axis0min,
00057                                    G4double          axis1min,
00058                                    G4double          axis0max,
00059                                    G4double          axis1max )
00060   : G4VTwistSurface(name, rot, tlate, 0, axis0, axis1,
00061                axis0min, axis1min, axis0max, axis1max)
00062 {   
00063    if (axis0 == kPhi && axis1 == kRho) {
00064       G4Exception("G4TwistTubsFlatSide::G4TwistTubsFlatSide()",
00065                   "GeomSolids0002", FatalErrorInArgument,
00066                   "Should swap axis0 and axis1!");
00067    }
00068    
00069    G4ThreeVector normal = rot.inverse()*n;
00070    fCurrentNormal.normal = normal.unit();   // in local coordinate system
00071    fIsValidNorm = true;
00072 
00073    SetCorners();
00074    SetBoundaries();
00075 
00076    fSurfaceArea = 1 ;  // not yet implemented. This is NOT a problem for tracking
00077 
00078 }
00079 
00080 
00081 
00082 G4TwistTubsFlatSide::G4TwistTubsFlatSide( const G4String        &name,
00083                                     G4double         EndInnerRadius[2],
00084                                     G4double         EndOuterRadius[2],
00085                                     G4double         DPhi,
00086                                     G4double         EndPhi[2],
00087                                     G4double         EndZ[2], 
00088                                     G4int            handedness ) 
00089   : G4VTwistSurface(name)
00090 {
00091    fHandedness = handedness;   // +z = +ve, -z = -ve
00092    fAxis[0]    = kRho;         // in local coordinate system
00093    fAxis[1]    = kPhi;
00094    G4int i     = (handedness < 0 ? 0 : 1);
00095    fAxisMin[0] = EndInnerRadius[i];  // Inner-hype radius at z=0
00096    fAxisMax[0] = EndOuterRadius[i];  // Outer-hype radius at z=0
00097    fAxisMin[1] = -0.5*DPhi;
00098    fAxisMax[1] = -fAxisMin[1];
00099    fCurrentNormal.normal.set(0, 0, (fHandedness < 0 ? -1 : 1)); 
00100          // Unit vector, in local coordinate system
00101    fRot.rotateZ(EndPhi[i]);
00102    fTrans.set(0, 0, EndZ[i]);
00103    fIsValidNorm = true;
00104 
00105    SetCorners();
00106    SetBoundaries();
00107 
00108    fSurfaceArea =  0.5*DPhi * (EndOuterRadius[i]*EndOuterRadius[i]  
00109                                - EndInnerRadius[i]*EndInnerRadius[i] ) ; 
00110 
00111 }
00112 
00113 
00114 //=====================================================================
00115 //* Fake default constructor ------------------------------------------
00116 
00117 G4TwistTubsFlatSide::G4TwistTubsFlatSide( __void__& a )
00118   : G4VTwistSurface(a), fSurfaceArea(0.)
00119 {
00120 }
00121 
00122 
00123 //=====================================================================
00124 //* destructor --------------------------------------------------------
00125 
00126 G4TwistTubsFlatSide::~G4TwistTubsFlatSide()
00127 {
00128 }
00129 
00130 //=====================================================================
00131 //* GetNormal ---------------------------------------------------------
00132 
00133 G4ThreeVector G4TwistTubsFlatSide::GetNormal(const G4ThreeVector & /* xx */ , 
00134                                              G4bool isGlobal)
00135 {
00136    if (isGlobal) {
00137       return ComputeGlobalDirection(fCurrentNormal.normal);
00138    } else {
00139       return fCurrentNormal.normal;
00140    }
00141 }
00142 
00143 //=====================================================================
00144 //* DistanceToSurface(p, v) -------------------------------------------
00145 
00146 G4int G4TwistTubsFlatSide::DistanceToSurface(const G4ThreeVector &gp,
00147                                        const G4ThreeVector &gv,
00148                                              G4ThreeVector  gxx[],
00149                                              G4double       distance[],
00150                                              G4int          areacode[],
00151                                              G4bool         isvalid[],
00152                                              EValidate      validate) 
00153 {
00154    fCurStatWithV.ResetfDone(validate, &gp, &gv);
00155    
00156    if (fCurStatWithV.IsDone()) {
00157       G4int i;
00158       for (i=0; i<fCurStatWithV.GetNXX(); i++) {
00159          gxx[i] = fCurStatWithV.GetXX(i);
00160          distance[i] = fCurStatWithV.GetDistance(i);
00161          areacode[i] = fCurStatWithV.GetAreacode(i);
00162          isvalid[i]  = fCurStatWithV.IsValid(i);
00163       }
00164       return fCurStatWithV.GetNXX();
00165    } else {
00166       // initialize
00167       G4int i;
00168       for (i=0; i<2; i++) {
00169          distance[i] = kInfinity;
00170          areacode[i] = sOutside;
00171          isvalid[i]  = false;
00172          gxx[i].set(kInfinity, kInfinity, kInfinity);
00173       }
00174    }
00175 
00176    G4ThreeVector p = ComputeLocalPoint(gp);
00177    G4ThreeVector v = ComputeLocalDirection(gv);
00178 
00179    //
00180    // special case!
00181    // if p is on surface, distance = 0. 
00182    //
00183 
00184    if (std::fabs(p.z()) == 0.) {   // if p is on the plane
00185       distance[0] = 0;
00186       G4ThreeVector xx = p;
00187       gxx[0] = ComputeGlobalPoint(xx);
00188       
00189       if (validate == kValidateWithTol) {
00190          areacode[0] = GetAreaCode(xx);
00191          if (!IsOutside(areacode[0])) {
00192             isvalid[0] = true;
00193          }
00194       } else if (validate == kValidateWithoutTol) {
00195          areacode[0] = GetAreaCode(xx, false);
00196          if (IsInside(areacode[0])) {
00197             isvalid[0] = true;
00198          }
00199       } else { // kDontValidate
00200          areacode[0] = sInside;
00201          isvalid[0] = true;
00202       }
00203 
00204       return 1;
00205    }
00206    //
00207    // special case end
00208    //
00209    
00210    if (v.z() == 0) { 
00211 
00212       fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0], 
00213                                      isvalid[0], 0, validate, &gp, &gv);
00214       return 0;
00215    }
00216    
00217    distance[0] = - (p.z() / v.z());
00218       
00219    G4ThreeVector xx = p + distance[0]*v;
00220    gxx[0] = ComputeGlobalPoint(xx);
00221 
00222    if (validate == kValidateWithTol) {
00223       areacode[0] = GetAreaCode(xx);
00224       if (!IsOutside(areacode[0])) {
00225          if (distance[0] >= 0) isvalid[0] = true;
00226       }
00227    } else if (validate == kValidateWithoutTol) {
00228       areacode[0] = GetAreaCode(xx, false);
00229       if (IsInside(areacode[0])) {
00230          if (distance[0] >= 0) isvalid[0] = true;
00231       }
00232    } else { // kDontValidate
00233       areacode[0] = sInside;
00234          if (distance[0] >= 0) isvalid[0] = true;
00235    }
00236 
00237    fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
00238                                   isvalid[0], 1, validate, &gp, &gv);
00239 
00240 #ifdef G4TWISTDEBUG
00241    G4cerr << "ERROR - G4TwistTubsFlatSide::DistanceToSurface(p,v)" << G4endl;
00242    G4cerr << "        Name        : " << GetName() << G4endl;
00243    G4cerr << "        xx          : " << xx << G4endl;
00244    G4cerr << "        gxx[0]      : " << gxx[0] << G4endl;
00245    G4cerr << "        dist[0]     : " << distance[0] << G4endl;
00246    G4cerr << "        areacode[0] : " << areacode[0] << G4endl;
00247    G4cerr << "        isvalid[0]  : " << isvalid[0]  << G4endl;
00248 #endif
00249    return 1;
00250 }
00251 
00252 //=====================================================================
00253 //* DistanceToSurface(p) ----------------------------------------------
00254 
00255 G4int G4TwistTubsFlatSide::DistanceToSurface(const G4ThreeVector &gp,
00256                                              G4ThreeVector  gxx[],
00257                                              G4double       distance[],
00258                                              G4int          areacode[])
00259 {
00260    // Calculate distance to plane in local coordinate,
00261    // then return distance and global intersection points.
00262    //  
00263 
00264    fCurStat.ResetfDone(kDontValidate, &gp);
00265 
00266    if (fCurStat.IsDone()) {
00267       G4int i;
00268       for (i=0; i<fCurStat.GetNXX(); i++) {
00269          gxx[i] = fCurStat.GetXX(i);
00270          distance[i] = fCurStat.GetDistance(i);
00271          areacode[i] = fCurStat.GetAreacode(i);
00272       }
00273       return fCurStat.GetNXX();
00274    } else {
00275       // initialize
00276       G4int i;
00277       for (i=0; i<2; i++) {
00278          distance[i] = kInfinity;
00279          areacode[i] = sOutside;
00280          gxx[i].set(kInfinity, kInfinity, kInfinity);
00281       }
00282    }
00283    
00284    G4ThreeVector p = ComputeLocalPoint(gp);
00285    G4ThreeVector xx;
00286 
00287    // The plane is placed on origin with making its normal 
00288    // parallel to z-axis. 
00289    if (std::fabs(p.z()) <= 0.5 * kCarTolerance) {   // if p is on the plane, return 1
00290       distance[0] = 0;
00291       xx = p;
00292    } else {
00293       distance[0] = std::fabs(p.z());
00294       xx.set(p.x(), p.y(), 0);  
00295    }
00296 
00297    gxx[0] = ComputeGlobalPoint(xx);
00298    areacode[0] = sInside;
00299    G4bool isvalid = true;
00300    fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
00301                              isvalid, 1, kDontValidate, &gp);
00302    return 1;
00303 
00304 }
00305 
00306 //=====================================================================
00307 //* GetAreaCode -------------------------------------------------------
00308 
00309 G4int G4TwistTubsFlatSide::GetAreaCode(const G4ThreeVector &xx, 
00310                                        G4bool withTol)
00311 {
00312 
00313    static const G4double rtol
00314      = 0.5*G4GeometryTolerance::GetInstance()->GetRadialTolerance();
00315    
00316    G4int areacode = sInside;
00317 
00318    if (fAxis[0] == kRho && fAxis[1] == kPhi) {
00319       G4int rhoaxis = 0;
00320      // G4int phiaxis = 0;
00321 
00322       G4ThreeVector dphimin;   // direction of phi-minimum boundary
00323       G4ThreeVector dphimax;   // direction of phi-maximum boundary
00324       dphimin = GetCorner(sC0Max1Min);
00325       dphimax = GetCorner(sC0Max1Max);   
00326       
00327       if (withTol) {
00328 
00329          G4bool isoutside = false;
00330 
00331          // test boundary of rho-axis
00332 
00333          if (xx.getRho() <= fAxisMin[rhoaxis] + rtol) {
00334 
00335             areacode |= (sAxis0 & (sAxisRho | sAxisMin)) | sBoundary; // rho-min
00336             if (xx.getRho() < fAxisMin[rhoaxis] - rtol) isoutside = true; 
00337             
00338          } else if (xx.getRho() >= fAxisMax[rhoaxis] - rtol) {
00339 
00340             areacode |= (sAxis0 & (sAxisRho | sAxisMax)) | sBoundary; // rho-max
00341             if (xx.getRho() > fAxisMax[rhoaxis] + rtol) isoutside = true; 
00342 
00343          }         
00344          
00345          // test boundary of phi-axis
00346 
00347          if (AmIOnLeftSide(xx, dphimin) >= 0) {           // xx is on dphimin
00348 
00349             areacode |= (sAxis1 & (sAxisPhi | sAxisMin)); 
00350             if   (areacode & sBoundary) areacode |= sCorner;  // xx is on the corner.
00351             else                        areacode |= sBoundary;
00352 
00353             if (AmIOnLeftSide(xx, dphimin) > 0) isoutside = true; 
00354 
00355          } else if (AmIOnLeftSide(xx, dphimax) <= 0) {    // xx is on dphimax
00356 
00357             areacode |= (sAxis1 & (sAxisPhi | sAxisMax)); 
00358             if   (areacode & sBoundary) areacode |= sCorner;  // xx is on the corner.
00359             else                        areacode |= sBoundary;
00360 
00361             if (AmIOnLeftSide(xx, dphimax) < 0) isoutside = true; 
00362 
00363          }
00364          
00365          // if isoutside = true, clear inside bit.
00366          // if not on boundary, add axis information. 
00367 
00368          if (isoutside) {
00369             G4int tmpareacode = areacode & (~sInside);
00370             areacode = tmpareacode;
00371          } else if ((areacode & sBoundary) != sBoundary) {
00372             areacode |= (sAxis0 & sAxisRho) | (sAxis1 & sAxisPhi);
00373          }
00374 
00375       } else {
00376 
00377          // out of boundary of rho-axis
00378 
00379          if (xx.getRho() < fAxisMin[rhoaxis]) {
00380             areacode |= (sAxis0 & (sAxisRho | sAxisMin)) | sBoundary;
00381          } else if (xx.getRho() > fAxisMax[rhoaxis]) {
00382             areacode |= (sAxis0 & (sAxisRho | sAxisMax)) | sBoundary;
00383          }
00384          
00385          // out of boundary of phi-axis
00386 
00387          if (AmIOnLeftSide(xx, dphimin, false) >= 0) {       // xx is leftside or
00388             areacode |= (sAxis1 & (sAxisPhi | sAxisMin)) ;   // boundary of dphimin
00389             if   (areacode & sBoundary) areacode |= sCorner; // xx is on the corner.
00390             else                        areacode |= sBoundary;
00391 
00392          } else if (AmIOnLeftSide(xx, dphimax, false) <= 0) { // xx is rightside or
00393             areacode |= (sAxis1 & (sAxisPhi | sAxisMax)) ;    // boundary of dphimax
00394             if   (areacode & sBoundary) areacode |= sCorner;  // xx is on the corner.
00395             else                        areacode |= sBoundary;
00396            
00397          }
00398 
00399          if ((areacode & sBoundary) != sBoundary) {
00400             areacode |= (sAxis0 & sAxisRho) | (sAxis1 & sAxisPhi);
00401          }
00402 
00403       }
00404       return areacode;
00405    } else {
00406 
00407       std::ostringstream message;
00408       message << "Feature NOT implemented !" << G4endl
00409               << "        fAxis[0] = " << fAxis[0] << G4endl
00410               << "        fAxis[1] = " << fAxis[1];
00411       G4Exception("G4TwistTubsFlatSide::GetAreaCode()", "GeomSolids0001",
00412                   FatalException, message);
00413    }
00414    return areacode;
00415 }
00416                                       
00417                                            
00418 //=====================================================================
00419 //* SetCorners --------------------------------------------------------
00420 
00421 void G4TwistTubsFlatSide::SetCorners()
00422 {
00423    // Set Corner points in local coodinate.
00424    
00425    if (fAxis[0] == kRho && fAxis[1] == kPhi) {
00426        
00427       G4int rhoaxis = 0;  // kRho
00428       G4int phiaxis = 1;  // kPhi
00429       
00430       G4double x, y, z;
00431       // corner of Axis0min and Axis1min
00432          x = fAxisMin[rhoaxis]*std::cos(fAxisMin[phiaxis]);
00433          y = fAxisMin[rhoaxis]*std::sin(fAxisMin[phiaxis]);
00434          z = 0;
00435          SetCorner(sC0Min1Min, x, y, z);
00436       // corner of Axis0max and Axis1min
00437          x = fAxisMax[rhoaxis]*std::cos(fAxisMin[phiaxis]);
00438          y = fAxisMax[rhoaxis]*std::sin(fAxisMin[phiaxis]);
00439          z = 0;
00440          SetCorner(sC0Max1Min, x, y, z);
00441       // corner of Axis0max and Axis1max
00442          x = fAxisMax[rhoaxis]*std::cos(fAxisMax[phiaxis]);
00443          y = fAxisMax[rhoaxis]*std::sin(fAxisMax[phiaxis]);
00444          z = 0;
00445          SetCorner(sC0Max1Max, x, y, z);
00446       // corner of Axis0min and Axis1max
00447          x = fAxisMin[rhoaxis]*std::cos(fAxisMax[phiaxis]);
00448          y = fAxisMin[rhoaxis]*std::sin(fAxisMax[phiaxis]);
00449          z = 0;
00450          SetCorner(sC0Min1Max, x, y, z);
00451        
00452    } else {
00453       std::ostringstream message;
00454       message << "Feature NOT implemented !" << G4endl
00455               << "        fAxis[0] = " << fAxis[0] << G4endl
00456               << "        fAxis[1] = " << fAxis[1];
00457       G4Exception("G4TwistTubsFlatSide::SetCorners()", "GeomSolids0001",
00458                   FatalException, message);
00459    }
00460 }
00461 
00462 //=====================================================================
00463 //* SetBoundaries() ---------------------------------------------------
00464 
00465 void G4TwistTubsFlatSide::SetBoundaries()
00466 {
00467    // Set direction-unit vector of phi-boundary-lines in local coodinate.
00468    // Don't call the function twice.
00469    
00470    if (fAxis[0] == kRho && fAxis[1] == kPhi) {
00471    
00472       G4ThreeVector direction;
00473       // sAxis0 & sAxisMin
00474       direction = GetCorner(sC0Min1Max) - GetCorner(sC0Min1Min);
00475       direction = direction.unit();
00476       SetBoundary(sAxis0 & (sAxisPhi | sAxisMin), direction,
00477                   GetCorner(sC0Min1Min), sAxisPhi);
00478                   
00479       // sAxis0 & sAxisMax
00480       direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min);
00481       direction = direction.unit();
00482       SetBoundary(sAxis0 & (sAxisPhi | sAxisMax), direction,
00483                   GetCorner(sC0Max1Min), sAxisPhi);
00484 
00485       // sAxis1 & sAxisMin
00486       direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
00487       direction = direction.unit();
00488       SetBoundary(sAxis1 & (sAxisRho | sAxisMin), direction,
00489                   GetCorner(sC0Min1Min), sAxisRho);
00490       
00491       // sAxis1 & sAxisMax
00492       direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max);
00493       direction = direction.unit();
00494       SetBoundary(sAxis1 & (sAxisRho | sAxisMax), direction,
00495                   GetCorner(sC0Min1Max), sAxisPhi);
00496    } else {
00497       std::ostringstream message;
00498       message << "Feature NOT implemented !" << G4endl
00499               << "        fAxis[0] = " << fAxis[0] << G4endl
00500               << "        fAxis[1] = " << fAxis[1];
00501       G4Exception("G4TwistTubsFlatSide::SetBoundaries()", "GeomSolids0001",
00502                   FatalException, message);
00503    }
00504 }
00505 
00506 //=====================================================================
00507 //* GetFacets() -------------------------------------------------------
00508 
00509 void G4TwistTubsFlatSide::GetFacets( G4int k, G4int n, G4double xyz[][3],
00510                                      G4int faces[][4], G4int iside ) 
00511 {
00512 
00513   G4ThreeVector p ;
00514 
00515   G4double rmin = fAxisMin[0] ;
00516   G4double rmax = fAxisMax[0] ;
00517   G4double phimin, phimax ;
00518 
00519   G4double r,phi ;
00520 
00521   G4int i,j ;
00522 
00523   G4int nnode,nface ;
00524 
00525   for ( i = 0 ; i<n ; i++ ) {
00526 
00527     r = rmin + i*(rmax-rmin)/(n-1) ;
00528 
00529     phimin = GetBoundaryMin(r) ; 
00530     phimax = GetBoundaryMax(r) ;
00531 
00532     for ( j = 0 ; j<k ; j++ )
00533     {
00534       phi = phimin + j*(phimax-phimin)/(k-1) ;
00535 
00536       nnode = GetNode(i,j,k,n,iside) ;
00537       p = SurfacePoint(phi,r,true) ;  // surface point in global coord.system
00538 
00539       xyz[nnode][0] = p.x() ;
00540       xyz[nnode][1] = p.y() ;
00541       xyz[nnode][2] = p.z() ;
00542 
00543       if ( i<n-1 && j<k-1 ) {   // conterclock wise filling
00544         
00545         nface = GetFace(i,j,k,n,iside) ;
00546 
00547         if (fHandedness < 0) {  // lower side 
00548           faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,-1) * ( GetNode(i  ,j  ,k,n,iside)+1) ;  
00549           faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,-1) * ( GetNode(i  ,j+1,k,n,iside)+1) ;
00550           faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,-1) * ( GetNode(i+1,j+1,k,n,iside)+1) ;
00551           faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,-1) * ( GetNode(i+1,j  ,k,n,iside)+1) ;
00552         } else {                // upper side
00553           faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,1) * ( GetNode(i  ,j  ,k,n,iside)+1) ;  
00554           faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,1) * ( GetNode(i+1,j  ,k,n,iside)+1) ;
00555           faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,1) * ( GetNode(i+1,j+1,k,n,iside)+1) ;
00556           faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,1) * ( GetNode(i  ,j+1,k,n,iside)+1) ;
00557 
00558         }
00559 
00560 
00561 
00562       }
00563     }
00564   }
00565 }

Generated on Mon May 27 17:50:04 2013 for Geant4 by  doxygen 1.4.7