G4ReplicaNavigation.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$
00028 //
00029 //
00030 // class G4ReplicaNavigation Implementation
00031 //
00032 // Author: P.Kent, 1996
00033 //
00034 // --------------------------------------------------------------------
00035 
00036 #include "G4ReplicaNavigation.hh"
00037 
00038 #include "G4AffineTransform.hh"
00039 #include "G4SmartVoxelProxy.hh"
00040 #include "G4SmartVoxelNode.hh"
00041 #include "G4VSolid.hh"
00042 #include "G4GeometryTolerance.hh"
00043 
00044 // ********************************************************************
00045 // Constructor
00046 // ********************************************************************
00047 //
00048 G4ReplicaNavigation::G4ReplicaNavigation()
00049 : fCheck(false), fVerbose(0)
00050 {
00051   kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
00052   kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
00053   kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
00054 }
00055 
00056 // ********************************************************************
00057 // Destructor
00058 // ********************************************************************
00059 //
00060 G4ReplicaNavigation::~G4ReplicaNavigation()
00061 {
00062 }
00063 
00064 // ********************************************************************
00065 // Inside
00066 // ********************************************************************
00067 //
00068 EInside
00069 G4ReplicaNavigation::Inside(const G4VPhysicalVolume *pVol,
00070                             const G4int replicaNo,
00071                             const G4ThreeVector &localPoint) const
00072 {
00073   EInside in = kOutside;
00074 
00075   // Replication data
00076   //
00077   EAxis axis;
00078   G4int nReplicas;
00079   G4double width, offset;
00080   G4bool consuming;
00081   
00082   G4double coord, rad2, rmin, tolRMax2, rmax, tolRMin2;
00083 
00084   pVol->GetReplicationData(axis, nReplicas, width, offset, consuming);
00085 
00086   switch (axis)
00087   {
00088     case kXAxis:
00089     case kYAxis:
00090     case kZAxis:
00091       coord = std::fabs(localPoint(axis))-width*0.5;
00092       if ( coord<=-kCarTolerance*0.5 )
00093       {
00094         in = kInside;
00095       }
00096       else if ( coord<=kCarTolerance*0.5 )
00097       {
00098         in = kSurface;
00099       }
00100       break;
00101     case kPhi:
00102       if ( localPoint.y()||localPoint.x() )
00103       {
00104         coord = std::fabs(std::atan2(localPoint.y(),localPoint.x()))-width*0.5;
00105         if ( coord<=-kAngTolerance*0.5 )
00106         {
00107           in = kInside;
00108         }
00109         else if ( coord<=kAngTolerance*0.5 )
00110         {
00111           in = kSurface;
00112         }
00113       }
00114       else
00115       {
00116         in = kSurface;
00117       }
00118       break;
00119     case kRho:
00120       rad2 = localPoint.perp2();
00121       rmax = (replicaNo+1)*width+offset;
00122       tolRMax2  = rmax-kRadTolerance*0.5;
00123       tolRMax2 *= tolRMax2;
00124       if ( rad2>tolRMax2 )
00125       {
00126         tolRMax2 = rmax+kRadTolerance*0.5;
00127         tolRMax2 *= tolRMax2;
00128         if ( rad2<=tolRMax2 )
00129         {
00130           in = kSurface;
00131         }
00132       }
00133       else
00134       {
00135         // Known to be inside outer radius
00136         //
00137         if ( replicaNo||offset )
00138         {
00139           rmin = rmax-width;
00140           tolRMin2 = rmin-kRadTolerance*0.5;
00141           tolRMin2 *= tolRMin2;
00142           if ( rad2>tolRMin2 )
00143           {
00144             tolRMin2 = rmin+kRadTolerance*0.5;
00145             tolRMin2 *= tolRMin2;
00146             if ( rad2>=tolRMin2 )
00147             {
00148               in = kInside;
00149             }
00150             else
00151             {
00152               in = kSurface;
00153             }
00154           }
00155         }
00156         else
00157         {
00158           in = kInside;
00159         }
00160       }
00161       break;
00162     default:
00163       G4Exception("G4ReplicaNavigation::Inside()", "GeomNav0002",
00164                   FatalException, "Unknown axis!");
00165       break;
00166   }
00167   return in;
00168 }
00169 
00170 // ********************************************************************
00171 // DistanceToOut
00172 // ********************************************************************
00173 //
00174 G4double
00175 G4ReplicaNavigation::DistanceToOut(const G4VPhysicalVolume *pVol,
00176                                    const G4int replicaNo,
00177                                    const G4ThreeVector &localPoint) const
00178 {
00179   // Replication data
00180   //
00181   EAxis axis;
00182   G4int nReplicas;
00183   G4double width,offset;
00184   G4bool consuming;
00185   
00186   G4double safety=0.;
00187   G4double safe1,safe2;
00188   G4double coord, rho, rmin, rmax;
00189 
00190   pVol->GetReplicationData(axis, nReplicas, width, offset, consuming);
00191   switch(axis)
00192   {
00193     case kXAxis:
00194     case kYAxis:
00195     case kZAxis:
00196        coord = localPoint(axis);
00197        safe1 = width*0.5-coord;
00198        safe2 = width*0.5+coord;
00199        safety = (safe1<=safe2) ? safe1 : safe2;
00200        break;
00201     case kPhi:
00202       if ( localPoint.y()<=0 )
00203       {
00204         safety = localPoint.x()*std::sin(width*0.5)
00205                + localPoint.y()*std::cos(width*0.5);
00206       }
00207       else
00208       {
00209         safety = localPoint.x()*std::sin(width*0.5)
00210                - localPoint.y()*std::cos(width*0.5);
00211       }
00212       break;
00213     case kRho:
00214       rho = localPoint.perp();
00215       rmax = width*(replicaNo+1)+offset;
00216       if ( replicaNo||offset )
00217       {
00218         rmin  = rmax-width;
00219         safe1 = rho-rmin;
00220         safe2 = rmax-rho;
00221         safety = (safe1<=safe2) ? safe1 : safe2;
00222       }
00223       else
00224       {
00225         safety = rmax-rho;
00226       }
00227       break;
00228     default:
00229      G4Exception("G4ReplicaNavigation::DistanceToOut()", "GeomNav0002",
00230                  FatalException, "Unknown axis!");
00231      break;
00232   }
00233   return (safety >= kCarTolerance) ? safety : 0;
00234 }
00235 
00236 // ********************************************************************
00237 // DistanceToOut
00238 // ********************************************************************
00239 //
00240 G4double
00241 G4ReplicaNavigation::DistanceToOut(const G4VPhysicalVolume *pVol,
00242                                    const G4int replicaNo,
00243                                    const G4ThreeVector &localPoint,
00244                                    const G4ThreeVector &localDirection,
00245                                    G4ExitNormal& arExitNormal ) const
00246 {
00247   // Replication data
00248   //
00249   EAxis axis;
00250   G4int nReplicas;
00251   G4double width, offset;
00252   G4bool consuming;
00253 
00254   G4double Dist=kInfinity;
00255   G4double coord, Comp, lindist;
00256   G4double signC = 0.0;
00257   G4ExitNormal candidateNormal; 
00258    
00259   static const G4ThreeVector VecCartAxes[3]=
00260    { G4ThreeVector(1.,0.,0.),G4ThreeVector(0.,1.,0.),G4ThreeVector(0.,0.,1.) };
00261   static G4ExitNormal::ESide SideCartAxesPlus [3]=
00262    { G4ExitNormal::kPX, G4ExitNormal::kPY, G4ExitNormal::kPZ };
00263   static G4ExitNormal::ESide SideCartAxesMinus[3]=
00264    { G4ExitNormal::kMX, G4ExitNormal::kMY, G4ExitNormal::kMZ };
00265 
00266   pVol->GetReplicationData(axis, nReplicas, width, offset, consuming);
00267   switch(axis)
00268   {
00269     case kXAxis:
00270     case kYAxis:
00271     case kZAxis:
00272       coord = localPoint(axis);
00273       Comp = localDirection(axis);
00274       if ( Comp>0 )
00275       {
00276         lindist = width*0.5-coord;
00277         Dist = (lindist>0) ? lindist/Comp : 0;
00278         signC= 1.0;
00279       }
00280       else if ( Comp<0 )
00281       {
00282         lindist = width*0.5+coord;
00283         Dist = (lindist>0) ? -lindist/Comp : 0;
00284         signC= -1.0;
00285       }
00286       else
00287       {
00288         Dist = kInfinity;
00289       }
00290       // signC = sign<G4double>(Comp)
00291       candidateNormal.exitNormal = ( signC * VecCartAxes[axis]);
00292       candidateNormal.calculated = true;
00293       candidateNormal.validConvex = true;
00294       candidateNormal.exitSide =
00295         (Comp>0) ? SideCartAxesPlus[axis] : SideCartAxesMinus[axis];
00296       break;
00297     case kPhi:
00298       Dist = DistanceToOutPhi(localPoint,localDirection,width,candidateNormal);
00299         // candidateNormal set in call
00300       break;
00301     case kRho:
00302       Dist = DistanceToOutRad(localPoint,localDirection,width,offset,
00303                               replicaNo,candidateNormal);
00304         // candidateNormal set in call
00305       break;
00306     default:
00307      G4Exception("G4ReplicaNavigation::DistanceToOut()", "GeomNav0002",
00308                  FatalException, "Unknown axis!");
00309      break;
00310   }
00311 
00312   arExitNormal= candidateNormal; // .exitNormal;
00313 
00314   return Dist;
00315 }
00316 
00317 // ********************************************************************
00318 // DistanceToOutPhi
00319 // ********************************************************************
00320 //
00321 G4double
00322 G4ReplicaNavigation::DistanceToOutPhi(const G4ThreeVector &localPoint,
00323                                       const G4ThreeVector &localDirection,
00324                                       const G4double width,
00325                                       G4ExitNormal& foundNormal ) const
00326 {
00327   // Phi Intersection
00328   // NOTE: width<=pi by definition
00329   //
00330   G4double sinSPhi= -2.0, cosSPhi= -2.0;
00331   G4double pDistS, pDistE, compS, compE, Dist, dist2, yi;
00332   G4ExitNormal::ESide sidePhi= G4ExitNormal::kNull;
00333   G4ThreeVector  candidateNormal;
00334 
00335   if ( (localPoint.x()!=0.0) || (localPoint.y()!=0.0) )
00336   {
00337     sinSPhi = std::sin(-width*0.5);  // SIN of starting phi plane
00338     cosSPhi = std::cos(width*0.5);   // COS of starting phi plane
00339 
00340     // pDist -ve when inside
00341     //
00342     pDistS = localPoint.x()*sinSPhi-localPoint.y()*cosSPhi;
00343      // Start plane at phi= -S
00344     pDistE = localPoint.x()*sinSPhi+localPoint.y()*cosSPhi;
00345      // End   plane at phi= +S
00346 
00347     // Comp -ve when in direction of outwards normal
00348     //
00349     compS = -sinSPhi*localDirection.x()+cosSPhi*localDirection.y();
00350     compE = -sinSPhi*localDirection.x()-cosSPhi*localDirection.y();
00351 
00352     if ( (pDistS<=0)&&(pDistE<=0) )
00353     {
00354       // Inside both phi *full* planes
00355       //
00356       if ( compS<0 )
00357       {
00358         dist2 = pDistS/compS;
00359         yi = localPoint.y()+dist2*localDirection.y();
00360 
00361         // Check intersecting with correct half-plane (no -> no intersect)
00362         //
00363         if ( yi<=0 )
00364         {
00365           Dist = (pDistS<=-kCarTolerance*0.5) ? dist2 : 0;
00366           sidePhi= G4ExitNormal::kSPhi; // tbc
00367         }
00368         else
00369         {
00370           Dist = kInfinity;
00371         }
00372       }
00373       else
00374       {
00375         Dist = kInfinity;
00376       }
00377       if ( compE<0 )
00378       {
00379         dist2 = pDistE/compE;
00380         
00381         // Only check further if < starting phi intersection
00382         //
00383         if ( dist2<Dist )
00384         {
00385           yi = localPoint.y()+dist2*localDirection.y();
00386 
00387           // Check intersecting with correct half-plane
00388           //
00389           if ( yi>=0 )
00390           {
00391             // Leaving via ending phi
00392             //
00393             Dist = (pDistE<=-kCarTolerance*0.5) ? dist2 : 0;
00394             sidePhi = G4ExitNormal::kEPhi;
00395           }
00396         }
00397       }
00398     }
00399     else if ( (pDistS>=0)&&(pDistE>=0) )
00400     {
00401       // Outside both *full* phi planes
00402       // if towards both >=0 then once inside will remain inside
00403       //
00404       Dist = ((compS>=0)&&(compE>=0)) ? kInfinity : 0;
00405     }
00406     else if ( (pDistS>0)&&(pDistE<0) )
00407     {
00408       // Outside full starting plane, inside full ending plane
00409       //
00410       if ( compE<0 )
00411       {      
00412         dist2 = pDistE/compE;
00413         yi = localPoint.y()+dist2*localDirection.y();
00414 
00415         // Check intersection in correct half-plane
00416         // (if not -> remain in extent)
00417         //
00418         Dist = (yi>0) ? dist2 : kInfinity;
00419         if( yi> 0 ) { sidePhi = G4ExitNormal::kEPhi; }
00420       }
00421       else  // Leaving immediately by starting phi
00422       {
00423         Dist = kInfinity;
00424       }
00425     }
00426     else
00427     {
00428       // Must be (pDistS<0)&&(pDistE>0)
00429       // Inside full starting plane, outside full ending plane
00430       //
00431       if ( compE>=0 )
00432       {
00433         if ( compS<0 )
00434         {
00435           dist2 = pDistS/compS;
00436           yi = localPoint.y()+dist2*localDirection.y();
00437 
00438           // Check intersection in correct half-plane
00439           // (if not -> remain in extent)
00440           //
00441           Dist = (yi<0) ? dist2 : kInfinity;
00442           if(yi<0)  { sidePhi = G4ExitNormal::kSPhi; }
00443         }
00444         else
00445         {
00446           Dist = kInfinity;
00447         }
00448       }
00449       else
00450       {
00451         // Leaving immediately by ending phi
00452         //
00453         Dist = 0;
00454         sidePhi= G4ExitNormal::kEPhi;
00455       }
00456     }
00457   }
00458   else
00459   {
00460     // On z axis + travel not || to z axis -> use direction vector
00461     //
00462     if( (std::fabs(localDirection.phi())<=width*0.5) )
00463     {
00464        Dist= kInfinity;
00465     }
00466     else
00467     {
00468        Dist= 0;
00469        sidePhi= G4ExitNormal::kMY;
00470     }
00471   }
00472 
00473   if(sidePhi == G4ExitNormal::kSPhi )
00474   {
00475     candidateNormal = G4ThreeVector(sinSPhi,-cosSPhi,0.) ;
00476   }
00477   else if (sidePhi == G4ExitNormal::kEPhi)
00478   {
00479     candidateNormal = G4ThreeVector(sinSPhi,cosSPhi,0.) ;
00480   }
00481   else if (sidePhi == G4ExitNormal::kMY )
00482   {
00483     candidateNormal = G4ThreeVector(0., -1.0, 0.); // Split -S and +S 'phi'
00484   }
00485   foundNormal.calculated= (sidePhi != G4ExitNormal::kNull );
00486   foundNormal.exitNormal= candidateNormal;
00487    
00488   return Dist;
00489 }
00490 
00491 // ********************************************************************
00492 // DistanceToOutRad
00493 // ********************************************************************
00494 //
00495 G4double
00496 G4ReplicaNavigation::DistanceToOutRad(const G4ThreeVector &localPoint,
00497                                       const G4ThreeVector &localDirection,
00498                                       const G4double width,
00499                                       const G4double offset,
00500                                       const G4int replicaNo,
00501                                       G4ExitNormal& foundNormal ) const
00502 {
00503   G4double rmin, rmax, t1, t2, t3, deltaR;
00504   G4double b, c, d2, srd;
00505   G4ExitNormal::ESide  sideR= G4ExitNormal::kNull;
00506    
00507   //
00508   // Radial Intersections
00509   //
00510   
00511   // Find intersction with cylinders at rmax/rmin
00512   // Intersection point (xi,yi,zi) on line
00513   // x=localPoint.x+t*localDirection.x etc.
00514   //
00515   // Intersects with x^2+y^2=R^2
00516   //
00517   // Hence (localDirection.x^2+localDirection.y^2)t^2+
00518   //     2t(localPoint.x*localDirection.x+localPoint.y*localDirection.y)+
00519   //        localPoint.x^2+localPoint.y^2-R^2=0
00520   //
00521   //            t1                t2                    t3
00522 
00523   rmin = replicaNo*width+offset;
00524   rmax = (replicaNo+1)*width+offset;
00525 
00526   t1 = 1.0-localDirection.z()*localDirection.z();   // since v normalised
00527   t2 = localPoint.x()*localDirection.x()+localPoint.y()*localDirection.y();
00528   t3 = localPoint.x()*localPoint.x()+localPoint.y()*localPoint.y();
00529   
00530   if ( t1>0 )        // Check not parallel
00531   {
00532     // Calculate srd, r exit distance
00533     //
00534     if ( t2>=0 )
00535     {
00536       // Delta r not negative => leaving via rmax
00537       //
00538       deltaR = t3-rmax*rmax;
00539     
00540       // NOTE: Should use
00541       // rho-rmax<-kRadTolerance*0.5 - [no sqrts for efficiency]
00542       //
00543       if ( deltaR<-kRadTolerance*0.5 )
00544       {
00545         b  = t2/t1;
00546         c  = deltaR/t1;
00547         srd = -b+std::sqrt(b*b-c);
00548         sideR= G4ExitNormal::kRMax;
00549       }
00550       else
00551       {
00552         // On tolerant boundary & heading outwards (or locally
00553         // perpendicular to) outer radial surface -> leaving immediately
00554         //
00555         srd = 0;
00556         sideR= G4ExitNormal::kRMax;
00557       }
00558     }
00559     else
00560     {
00561       // Possible rmin intersection
00562       //
00563       if (rmin)
00564       {
00565         deltaR = t3-rmin*rmin;
00566         b  = t2/t1;
00567         c  = deltaR/t1;
00568         d2 = b*b-c;
00569         if ( d2>=0 )
00570         {
00571           // Leaving via rmin
00572           // NOTE: Should use
00573           // rho-rmin>kRadTolerance*0.5 - [no sqrts for efficiency]
00574           //
00575           srd = (deltaR>kRadTolerance*0.5) ? -b-std::sqrt(d2) : 0;
00576           // Is the following more accurate ?   - called 'issue' below
00577           // srd = (deltaR>kRadTolerance*0.5) ? c/( -b - std::sqrt(d2)) : 0.0;
00578           sideR= G4ExitNormal::kRMin;
00579         }
00580         else
00581         {
00582           // No rmin intersect -> must be rmax intersect
00583           //
00584           deltaR = t3-rmax*rmax;
00585           c  = deltaR/t1;
00586           srd = -b+std::sqrt(b*b-c); //  See issue above
00587           sideR= G4ExitNormal::kRMax;
00588         }
00589       }
00590       else
00591       {
00592         // No rmin intersect -> must be rmax intersect
00593         //
00594         deltaR = t3-rmax*rmax;
00595         b  = t2/t1;
00596         c  = deltaR/t1;
00597         srd = -b+std::sqrt(b*b-c);  // See issue above
00598         sideR= G4ExitNormal::kRMax;
00599       }
00600     }
00601   }
00602   else
00603   {
00604     srd=kInfinity;
00605     sideR= G4ExitNormal::kNull;
00606   }
00607    
00608   if( sideR != G4ExitNormal::kNull ) // if ((side == kRMax) || (side==kRMin))
00609   {
00610      // Note: returned vector not explicitly normalised
00611      // (divided by fRMax for unit vector)
00612      G4double xi, yi;
00613      xi = localPoint.x() + srd*localDirection.x() ;
00614      yi = localPoint.y() + srd*localDirection.y() ;
00615      G4ThreeVector normalR = G4ThreeVector(xi,yi,0.0);
00616      
00617      if( sideR == G4ExitNormal::kRMax ){
00618         normalR *=   1.0/rmax;
00619      }else{
00620         normalR *= (-1.0)/rmin;
00621      }
00622      foundNormal.exitNormal= normalR;
00623      foundNormal.calculated= true;
00624      foundNormal.validConvex = (sideR == G4ExitNormal::kRMax);
00625      foundNormal.exitSide = sideR;
00626   }else{
00627      foundNormal.calculated= false;
00628   }
00629    
00630   return srd;
00631 }
00632 
00633 // ********************************************************************
00634 // ComputeTransformation
00635 //
00636 // Setup transformation and transform point into local system
00637 // ********************************************************************
00638 //
00639 void
00640 G4ReplicaNavigation::ComputeTransformation(const G4int replicaNo,
00641                                                  G4VPhysicalVolume* pVol,
00642                                                  G4ThreeVector& point) const
00643 {
00644   G4double val,cosv,sinv,tmpx,tmpy;
00645 
00646   // Replication data
00647   //
00648   EAxis axis;
00649   G4int nReplicas;
00650   G4double width,offset;
00651   G4bool consuming;
00652 
00653   pVol->GetReplicationData(axis, nReplicas, width, offset, consuming);
00654 
00655   switch (axis)
00656   {
00657     case kXAxis:
00658       val = -width*0.5*(nReplicas-1)+width*replicaNo;
00659       pVol->SetTranslation(G4ThreeVector(val,0,0));
00660       point.setX(point.x()-val);
00661       break;
00662     case kYAxis:
00663       val = -width*0.5*(nReplicas-1)+width*replicaNo;
00664       pVol->SetTranslation(G4ThreeVector(0,val,0));
00665       point.setY(point.y()-val);
00666       break;
00667     case kZAxis:
00668       val = -width*0.5*(nReplicas-1)+width*replicaNo;
00669       pVol->SetTranslation(G4ThreeVector(0,0,val));
00670       point.setZ(point.z()-val);
00671       break;
00672     case kPhi:
00673       val = -(offset+width*(replicaNo+0.5));
00674       SetPhiTransformation(val,pVol);
00675       cosv = std::cos(val);
00676       sinv = std::sin(val);
00677       tmpx = point.x()*cosv-point.y()*sinv;
00678       tmpy = point.x()*sinv+point.y()*cosv;
00679       point.setY(tmpy);
00680       point.setX(tmpx);
00681       break;
00682     case kRho:
00683       // No setup required for radial case
00684     default:
00685       break;
00686   }
00687 }
00688 
00689 // ********************************************************************
00690 // ComputeTransformation
00691 //
00692 // Setup transformation into local system
00693 // ********************************************************************
00694 //
00695 void
00696 G4ReplicaNavigation::ComputeTransformation(const G4int replicaNo,
00697                                                  G4VPhysicalVolume* pVol) const
00698 {
00699   G4double val;
00700 
00701   // Replication data
00702   //
00703   EAxis axis;
00704   G4int nReplicas;
00705   G4double width, offset;
00706   G4bool consuming;
00707 
00708   pVol->GetReplicationData(axis, nReplicas, width, offset, consuming);
00709 
00710   switch (axis)
00711   {
00712     case kXAxis:
00713       val = -width*0.5*(nReplicas-1)+width*replicaNo;
00714       pVol->SetTranslation(G4ThreeVector(val,0,0));
00715       break;
00716     case kYAxis:
00717       val = -width*0.5*(nReplicas-1)+width*replicaNo;
00718       pVol->SetTranslation(G4ThreeVector(0,val,0));
00719       break;
00720     case kZAxis:
00721       val = -width*0.5*(nReplicas-1)+width*replicaNo;
00722       pVol->SetTranslation(G4ThreeVector(0,0,val));
00723       break;
00724     case kPhi:
00725       val = -(offset+width*(replicaNo+0.5));
00726       SetPhiTransformation(val,pVol);
00727       break;
00728     case kRho:
00729       // No setup required for radial case
00730     default:
00731       break;
00732   }
00733 }
00734 
00735 // ********************************************************************
00736 // ComputeStep
00737 // ********************************************************************
00738 //
00739 G4double
00740 G4ReplicaNavigation::ComputeStep(const G4ThreeVector &globalPoint,
00741                                  const G4ThreeVector &globalDirection,
00742                                  const G4ThreeVector &localPoint,
00743                                  const G4ThreeVector &localDirection,
00744                                  const G4double currentProposedStepLength,
00745                                        G4double &newSafety,
00746                                        G4NavigationHistory &history,
00747                                  // std::pair<G4bool,G4bool> &validAndCalculated
00748                                        G4bool &validExitNormal,
00749                                        G4bool &calculatedExitNormal, 
00750                                        G4ThreeVector &exitNormalVector,
00751                                        G4bool &exiting,
00752                                        G4bool &entering,
00753                                        G4VPhysicalVolume *(*pBlockedPhysical),
00754                                        G4int &blockedReplicaNo )
00755 {
00756   G4VPhysicalVolume *repPhysical, *motherPhysical;
00757   G4VPhysicalVolume *samplePhysical, *blockedExitedVol=0;
00758   G4LogicalVolume *repLogical;
00759   G4VSolid *motherSolid;
00760   G4ThreeVector repPoint, repDirection, sampleDirection;
00761   G4double ourStep=currentProposedStepLength;
00762   G4double ourSafety=kInfinity;
00763   G4double sampleStep, sampleSafety, motherStep, motherSafety;
00764   G4int localNoDaughters, sampleNo;
00765   G4int depth;
00766   G4ExitNormal exitNormalStc;
00767   // G4int depthDeterminingStep= -1; // Useful only for debugging - for now
00768 
00769   calculatedExitNormal= false;
00770   
00771   // Exiting normal optimisation
00772   //
00773   if ( exiting&&validExitNormal )
00774   {
00775     if ( localDirection.dot(exitNormalVector)>=kMinExitingNormalCosine )
00776     {
00777       // Block exited daughter volume
00778       //
00779       blockedExitedVol = *pBlockedPhysical;
00780       ourSafety = 0;
00781     }
00782   }
00783   exiting  = false;
00784   entering = false;
00785 
00786   repPhysical = history.GetTopVolume();
00787   repLogical  = repPhysical->GetLogicalVolume();
00788 
00789   //
00790   // Compute intersection with replica boundaries & replica safety
00791   //
00792 
00793   sampleSafety = DistanceToOut(repPhysical,
00794                                history.GetTopReplicaNo(),
00795                                localPoint);
00796   G4ExitNormal normalOutStc;
00797   const G4int topDepth= history.GetDepth();
00798 
00799   if ( sampleSafety<ourSafety )
00800   {
00801     ourSafety = sampleSafety;
00802   }
00803   if ( sampleSafety<ourStep )
00804   {
00805 
00806     sampleStep = DistanceToOut(repPhysical,
00807                                history.GetTopReplicaNo(),
00808                                localPoint,
00809                                localDirection,
00810                                normalOutStc);
00811     if ( sampleStep<ourStep )
00812     {
00813       ourStep = sampleStep;
00814       exiting = true;
00815       validExitNormal = normalOutStc.validConvex; // false; -> Old,Conservative
00816 
00817       exitNormalStc= normalOutStc;
00818       exitNormalStc.exitNormal= history.GetTopTransform().Inverse().
00819                                 TransformAxis(normalOutStc.exitNormal);
00820       calculatedExitNormal= true;
00821     }
00822   }
00823   const G4int secondDepth= topDepth;
00824   depth = secondDepth;
00825   
00826   while ( history.GetVolumeType(depth)==kReplica )
00827   {
00828     const G4AffineTransform& GlobalToLocal= history.GetTransform(depth);
00829     repPoint    = GlobalToLocal.TransformPoint(globalPoint);
00830     // repPoint = history.GetTransform(depth).TransformPoint(globalPoint);
00831  
00832     sampleSafety = DistanceToOut(history.GetVolume(depth),
00833                                  history.GetReplicaNo(depth),
00834                                  repPoint);
00835     if ( sampleSafety < ourSafety )
00836     {
00837       ourSafety = sampleSafety;
00838     }
00839     if ( sampleSafety < ourStep )
00840     {
00841       G4ThreeVector newLocalDirection = GlobalToLocal.TransformAxis(globalDirection);
00842       sampleStep = DistanceToOut(history.GetVolume(depth),
00843                                  history.GetReplicaNo(depth),
00844                                  repPoint,
00845                                  newLocalDirection,
00846                                  normalOutStc);
00847       if ( sampleStep < ourStep )
00848       {
00849         ourStep = sampleStep;
00850         exiting = true;
00851        
00852         // As step is limited by this level, must set Exit Normal
00853         //
00854         G4ThreeVector localExitNorm= normalOutStc.exitNormal;
00855         G4ThreeVector globalExitNorm=
00856             GlobalToLocal.Inverse().TransformAxis(localExitNorm);
00857 
00858         exitNormalStc= normalOutStc; // Normal, convex, calculated, side
00859         exitNormalStc.exitNormal= globalExitNorm;
00860         calculatedExitNormal= true;
00861       }
00862     }
00863     depth--;
00864   }
00865  
00866   // Compute mother safety & intersection
00867   //
00868   G4ThreeVector exitVectorMother;
00869   G4bool        exitConvex= false; // Value obtained in DistanceToOut(p,v) call
00870   G4ExitNormal  motherNormalStc;
00871 
00872   repPoint = history.GetTransform(depth).TransformPoint(globalPoint);
00873   motherPhysical = history.GetVolume(depth);
00874   motherSolid = motherPhysical->GetLogicalVolume()->GetSolid();
00875   motherSafety = motherSolid->DistanceToOut(repPoint);
00876   repDirection = history.GetTransform(depth).TransformAxis(globalDirection);
00877 
00878   motherStep = motherSolid->DistanceToOut(repPoint,repDirection,true,
00879                                           &exitConvex,&exitVectorMother);
00880   if( exitConvex )
00881   {
00882      motherNormalStc = G4ExitNormal( exitVectorMother, true, false,
00883                                      G4ExitNormal::kMother);
00884      calculatedExitNormal= true;
00885   }
00886   const G4AffineTransform& globalToLocalTop = history.GetTopTransform();
00887 
00888   G4bool motherDeterminedStep= (motherStep<ourStep);
00889 
00890   if( (!exitConvex) && motherDeterminedStep )
00891   {
00892      exitVectorMother= motherSolid->SurfaceNormal( repPoint );
00893      motherNormalStc= G4ExitNormal( exitVectorMother, true, false,
00894                                     G4ExitNormal::kMother);
00895      // CalculatedExitNormal -> true;
00896      // Convex               -> false: do not know value
00897      // ExitSide             -> kMother (or kNull)
00898  
00899      calculatedExitNormal= true;
00900   }
00901   if( motherDeterminedStep)
00902   {
00903      G4ThreeVector globalExitNormalTop=
00904         globalToLocalTop.Inverse().TransformAxis(exitVectorMother);
00905      
00906      exitNormalStc= motherNormalStc;
00907      exitNormalStc.exitNormal= globalExitNormalTop;
00908   }
00909 
00910   // Push in principle no longer necessary. G4Navigator now takes care of ...
00911   // Removing this will however generate warnings for pushed particles from
00912   // G4Navigator, particularly for the case of 3D replicas (Cartesian or
00913   // combined Radial/Phi cases).
00914   // Requires further investigation and eventually reimplementation of
00915   // LevelLocate() to take into account point and direction ...
00916   //
00917   if  ( ( !ourStep && (sampleSafety<0.5*kCarTolerance) )
00918      && ( repLogical->GetSolid()->Inside(localPoint)==kSurface ) )
00919   {
00920     ourStep += kCarTolerance;
00921   }
00922 
00923   if ( motherSafety<ourSafety )
00924   {
00925     ourSafety = motherSafety;
00926   }
00927 
00928 #ifdef G4VERBOSE
00929   if ( fCheck )
00930   {
00931     if( motherSolid->Inside(localPoint)==kOutside )
00932     {
00933       std::ostringstream message;
00934       message << "Point outside volume !" << G4endl
00935               << "          Point " << localPoint
00936               << " is outside current volume " << motherPhysical->GetName()
00937               << G4endl;
00938       G4double estDistToSolid= motherSolid->DistanceToIn(localPoint); 
00939       message << "          Estimated isotropic distance to solid (distToIn)= " 
00940               << estDistToSolid << G4endl;
00941       if( estDistToSolid > 100.0 * kCarTolerance )
00942       {
00943         motherSolid->DumpInfo();
00944         G4Exception("G4ReplicaNavigation::ComputeStep()",
00945                     "GeomNav0003", FatalException, message,
00946                     "Point is far outside Current Volume !" ); 
00947       }
00948       else
00949         G4Exception("G4ReplicaNavigation::ComputeStep()",
00950                     "GeomNav1002", JustWarning, message,
00951                     "Point is a little outside Current Volume."); 
00952     }
00953   }
00954 #endif
00955 
00956   // Comparison of steps may need precision protection
00957   //
00958 #if 1
00959   if( motherDeterminedStep)
00960   {
00961     ourStep = motherStep;
00962     exiting = true;
00963   }
00964 
00965   // Transform it to the Grand-Mother Reference Frame (current convention)
00966   //
00967   if ( calculatedExitNormal )
00968   {
00969     if ( motherDeterminedStep )
00970     {
00971       exitNormalVector= motherNormalStc.exitNormal;
00972     }else{
00973       G4ThreeVector exitNormalGlobal= exitNormalStc.exitNormal;
00974       exitNormalVector= globalToLocalTop.TransformAxis(exitNormalGlobal);
00975       // exitNormalVector= globalToLocal2nd.TransformAxis(exitNormalGlobal);
00976       // Alt Make it in one go to Grand-Mother, avoiding transform below
00977     }
00978     // Transform to Grand-mother reference frame
00979     const G4RotationMatrix* rot = motherPhysical->GetRotation();
00980     if ( rot )
00981     {
00982       exitNormalVector *= rot->inverse();
00983     }
00984 
00985   }
00986   else
00987   {
00988     validExitNormal = false;
00989   }
00990 
00991 #else
00992   if ( motherSafety<=ourStep )
00993   {
00994     if ( motherStep<=ourStep )
00995     {
00996       ourStep = motherStep;
00997       exiting = true;
00998       if ( validExitNormal )
00999       {
01000         const G4RotationMatrix* rot = motherPhysical->GetRotation();
01001         if ( rot )
01002         {
01003           exitNormal *= rot->inverse();
01004         }
01005       }
01006     }
01007     else
01008     {
01009       validExitNormal = false;
01010       // calculatedExitNormal= false;
01011     }
01012   }
01013 #endif
01014 
01015 
01016   G4bool daughterDeterminedStep=false;
01017   G4ThreeVector daughtNormRepCrd;
01018      // Exit normal of daughter transformed to
01019      // the coordinate system of Replica (i.e. last depth)
01020 
01021   //
01022   // Compute daughter safeties & intersections
01023   //
01024   localNoDaughters = repLogical->GetNoDaughters();
01025   for ( sampleNo=localNoDaughters-1; sampleNo>=0; sampleNo-- )
01026   {
01027     samplePhysical = repLogical->GetDaughter(sampleNo);
01028     if ( samplePhysical!=blockedExitedVol )
01029     {
01030       G4ThreeVector localExitNorm;
01031       G4ThreeVector normReplicaCoord;
01032 
01033       G4AffineTransform sampleTf(samplePhysical->GetRotation(),
01034                                  samplePhysical->GetTranslation());
01035       sampleTf.Invert();
01036       const G4ThreeVector samplePoint =
01037                         sampleTf.TransformPoint(localPoint);
01038       const G4VSolid* sampleSolid =
01039                         samplePhysical->GetLogicalVolume()->GetSolid();
01040       const G4double sampleSafetyDistance =
01041                         sampleSolid->DistanceToIn(samplePoint);
01042       if ( sampleSafetyDistance<ourSafety )
01043       {
01044         ourSafety = sampleSafetyDistance;
01045       }
01046       if ( sampleSafetyDistance<=ourStep )
01047       {
01048         sampleDirection = sampleTf.TransformAxis(localDirection);
01049         const G4double sampleStepDistance =
01050                         sampleSolid->DistanceToIn(samplePoint,sampleDirection);
01051         if ( sampleStepDistance<=ourStep )
01052         {
01053           daughterDeterminedStep= true;
01054 
01055           ourStep  = sampleStepDistance;
01056           entering = true;
01057           exiting  = false;
01058           *pBlockedPhysical = samplePhysical;
01059           blockedReplicaNo  = -1;
01060 
01061 #ifdef DAUGHTER_NORMAL_ALSO
01062           // This norm can be calculated later, if needed daughter is available
01063           localExitNorm = sampleSolid->SurfaceNormal(samplePoint);
01064           daughtNormRepCrd = sampleTf.Inverse().TransformAxis(localExitNorm);
01065 #endif
01066           
01067 #ifdef G4VERBOSE
01068           // Check to see that the resulting point is indeed in/on volume.
01069           // This check could eventually be made only for successful candidate.
01070 
01071           if ( ( fCheck ) && ( sampleStepDistance < kInfinity ) )
01072           {
01073             G4ThreeVector intersectionPoint;
01074             intersectionPoint= samplePoint
01075                              + sampleStepDistance * sampleDirection;
01076             EInside insideIntPt= sampleSolid->Inside(intersectionPoint); 
01077             if ( insideIntPt != kSurface )
01078             {
01079               G4int oldcoutPrec = G4cout.precision(16); 
01080               std::ostringstream message;
01081               message << "Navigator gets conflicting response from Solid."
01082                       << G4endl
01083                       << "          Inaccurate DistanceToIn for solid "
01084                       << sampleSolid->GetName() << G4endl
01085                       << "          Solid gave DistanceToIn = "
01086                       << sampleStepDistance << " yet returns " ;
01087               if ( insideIntPt == kInside )
01088                 message << "-kInside-"; 
01089               else if ( insideIntPt == kOutside )
01090                 message << "-kOutside-";
01091               else
01092                 message << "-kSurface-"; 
01093               message << " for this point !" << G4endl
01094                       << "          Point = " << intersectionPoint << G4endl;
01095               if ( insideIntPt != kInside )
01096                 message << "        DistanceToIn(p) = " 
01097                        << sampleSolid->DistanceToIn(intersectionPoint)
01098                        << G4endl;
01099               if ( insideIntPt != kOutside ) 
01100                 message << "        DistanceToOut(p) = " 
01101                        << sampleSolid->DistanceToOut(intersectionPoint);
01102               G4Exception("G4ReplicaNavigation::ComputeStep()", 
01103                           "GeomNav1002", JustWarning, message); 
01104               G4cout.precision(oldcoutPrec);
01105             }
01106           }
01107 #endif
01108         }
01109       }
01110     }
01111   }
01112 
01113   calculatedExitNormal &= (!daughterDeterminedStep);
01114 
01115 #ifdef DAUGHTER_NORMAL_ALSO
01116   if( daughterDeterminedStep )
01117   {
01118     // G4ThreeVector daughtNormGlobal =
01119     //   GlobalToLastDepth.Inverse().TransformAxis(daughtNormRepCrd);
01120     // ==> Can calculate it, but have no way to transmit it to caller (for now)
01121 
01122     exitNormalVector=globalToLocalTop.Inverse().TransformAxis(daughtNormGlobal);
01123     validExitNormal= false; // Entering daughter - never convex for parent
01124 
01125     calculatedExitNormal= true;
01126   }
01127   // calculatedExitNormal= true;  // Force it to true -- dubious
01128 #endif
01129 
01130   newSafety = ourSafety;
01131   return ourStep;
01132 }
01133 
01134 // ********************************************************************
01135 // ComputeSafety
01136 //
01137 // Compute the isotropic distance to current volume's boundaries
01138 // and to daughter volumes.
01139 // ********************************************************************
01140 //
01141 G4double
01142 G4ReplicaNavigation::ComputeSafety(const G4ThreeVector &globalPoint,
01143                                    const G4ThreeVector &localPoint,
01144                                          G4NavigationHistory &history,
01145                                    const G4double )
01146 {
01147   G4VPhysicalVolume *repPhysical, *motherPhysical;
01148   G4VPhysicalVolume *samplePhysical, *blockedExitedVol=0;
01149   G4LogicalVolume *repLogical;
01150   G4VSolid *motherSolid;
01151   G4ThreeVector repPoint;
01152   G4double ourSafety=kInfinity;
01153   G4double sampleSafety;
01154   G4int localNoDaughters, sampleNo;
01155   G4int depth;
01156 
01157   repPhysical = history.GetTopVolume();
01158   repLogical  = repPhysical->GetLogicalVolume();
01159 
01160   //
01161   // Compute intersection with replica boundaries & replica safety
01162   //
01163 
01164   sampleSafety = DistanceToOut(history.GetTopVolume(),
01165                                history.GetTopReplicaNo(),
01166                                localPoint);
01167   if ( sampleSafety<ourSafety )
01168   {
01169     ourSafety = sampleSafety;
01170   }
01171 
01172   depth = history.GetDepth()-1;
01173   while ( history.GetVolumeType(depth)==kReplica )
01174   {
01175     repPoint = history.GetTransform(depth).TransformPoint(globalPoint);
01176     sampleSafety = DistanceToOut(history.GetVolume(depth),
01177                                  history.GetReplicaNo(depth),
01178                                  repPoint);
01179     if ( sampleSafety<ourSafety )
01180     {
01181       ourSafety = sampleSafety;
01182     }
01183     depth--;
01184   }
01185 
01186   // Compute mother safety & intersection
01187   //
01188   repPoint = history.GetTransform(depth).TransformPoint(globalPoint);
01189   motherPhysical = history.GetVolume(depth);
01190   motherSolid = motherPhysical->GetLogicalVolume()->GetSolid();
01191   sampleSafety = motherSolid->DistanceToOut(repPoint);
01192 
01193   if ( sampleSafety<ourSafety )
01194   {
01195     ourSafety = sampleSafety;
01196   }
01197 
01198   // Compute daughter safeties & intersections
01199   //
01200   localNoDaughters = repLogical->GetNoDaughters();
01201   for ( sampleNo=localNoDaughters-1; sampleNo>=0; sampleNo-- )
01202   {
01203     samplePhysical = repLogical->GetDaughter(sampleNo);
01204     if ( samplePhysical!=blockedExitedVol )
01205     {
01206       G4AffineTransform sampleTf(samplePhysical->GetRotation(),
01207                                  samplePhysical->GetTranslation());
01208       sampleTf.Invert();
01209       const G4ThreeVector samplePoint =
01210                             sampleTf.TransformPoint(localPoint);
01211       const G4VSolid *sampleSolid =
01212                             samplePhysical->GetLogicalVolume()->GetSolid();
01213       const G4double sampleSafetyDistance =
01214                             sampleSolid->DistanceToIn(samplePoint);
01215       if ( sampleSafetyDistance<ourSafety )
01216       {
01217         ourSafety = sampleSafetyDistance;
01218       }
01219     }
01220   }
01221   return ourSafety;
01222 }
01223 
01224 // ********************************************************************
01225 // BackLocate
01226 // ********************************************************************
01227 //
01228 EInside
01229 G4ReplicaNavigation::BackLocate(G4NavigationHistory &history,
01230                           const G4ThreeVector &globalPoint,
01231                                 G4ThreeVector &localPoint,
01232                           const G4bool &exiting,
01233                                 G4bool &notKnownInside ) const
01234 {
01235   G4VPhysicalVolume *pNRMother=0;
01236   G4VSolid *motherSolid;
01237   G4ThreeVector repPoint, goodPoint;
01238   G4int mdepth, depth, cdepth;
01239   EInside insideCode;
01240 
01241   cdepth = history.GetDepth();
01242   
01243   // Find non replicated mother
01244   //
01245   for ( mdepth=cdepth-1; mdepth>=0; mdepth-- )
01246   {
01247     if ( history.GetVolumeType(mdepth)!=kReplica )
01248     {
01249       pNRMother = history.GetVolume(mdepth);
01250       break;
01251     }
01252   }
01253 
01254   if( pNRMother==0 ) 
01255   {
01256     // All the tree of mother volumes were Replicas. 
01257     // This is an error, as the World volume must be a Placement
01258     //
01259     G4Exception("G4ReplicaNavigation::BackLocate()", "GeomNav0002",
01260                 FatalException, "The World volume must be a Placement!");
01261     return kInside;
01262   }
01263 
01264   motherSolid = pNRMother->GetLogicalVolume()->GetSolid();
01265   goodPoint = history.GetTransform(mdepth).TransformPoint(globalPoint);
01266   insideCode = motherSolid->Inside(goodPoint);
01267   if ( (insideCode==kOutside)||((insideCode==kSurface)&&exiting) )
01268   {
01269     // Outside mother -> back up to mother level
01270     // Locate.. in Navigator will back up one more level
01271     // localPoint not required
01272     //
01273     history.BackLevel(cdepth-mdepth);
01274     //      localPoint = goodPoint;
01275   }
01276   else
01277   {
01278     notKnownInside = false;
01279 
01280     // Still within replications
01281     // Check down: if on outside stop at this level
01282     //
01283     for ( depth=mdepth+1; depth<cdepth; depth++)
01284     {
01285       repPoint = history.GetTransform(depth).TransformPoint(globalPoint);
01286       insideCode = Inside(history.GetVolume(depth),
01287                           history.GetReplicaNo(depth),
01288                           repPoint);
01289       if ( (insideCode==kOutside)||((insideCode==kSurface)&&exiting) )
01290       {
01291         localPoint = goodPoint;
01292         history.BackLevel(cdepth-depth);
01293         return insideCode;
01294       }
01295       else
01296       {
01297         goodPoint = repPoint;
01298       }
01299     }
01300     localPoint = history.GetTransform(depth).TransformPoint(globalPoint);
01301     insideCode = Inside(history.GetVolume(depth),
01302                         history.GetReplicaNo(depth),
01303                         localPoint);
01304     // If outside level, set localPoint = coordinates in reference system
01305     // of *previous* level - location code in navigator will back up one
01306     // level [And also manage blocking]
01307     //
01308     if ( (insideCode==kOutside)||((insideCode==kSurface)&&exiting) )
01309     {
01310       localPoint = goodPoint;
01311     }
01312   }
01313   return insideCode;
01314 }

Generated on Mon May 27 17:49:43 2013 for Geant4 by  doxygen 1.4.7