Flugg/WrapNorml.cxx

   1
   2 // Flugg tag
   3
   4 ////////////////////////////////////////////////////////////////////
   5 //
   6 // WrapNorml.hh - Sara Vanini
   7 //
   8 // Wrapper for computing normal unit-vector in global coordinates.
   9 //
  10 // Fluka requires normal vector exiting from final position (of
  11 // particle) volume, that is: entering in volume of initial position.
  12 // Geant4 always computes the normal vector exiting from the volume.
  13 // In GetLocalExitNormal() call the volume is the pre-step volume
  14 // (so G4 normal vector sign is opposite of fluka-required normal).
  15 // If IsLocalExitNormalValid=false, normal value is computed from
  16 // init-step volume (in this case sign must change), or from
  17 // end-step volume (the sign is the same). Normal vector is computed
  18 // always on boundary between volumes, in global coordinates (to take
  19 // rotation of parameterised volumes in hierarchy in consideration).
  20 // So: nrmlwr returns inwards pointing unit normal of the shape for
  21 // surface closest to the point returned by the navigator (last step
  22 // end-point).
  23 //
  24 // modified 10/III/99
  25 // modified 25/V/00
  26 // modified 7/VI/00 for boundary-crossing in case of relocation
  27 // modified 5/VII/00 geometry error on boundary fixed
  28 // modified 24.10.01: by I. Hrivnacova
  29 //   functions declarations separated from implementation
  30 //   (moved to Wrappers.hh);
  31 //
  32 ////////////////////////////////////////////////////////////////////
  33
  34 #include "Wrappers.hh"
  35 #include "WrapUtils.hh"
  36 #include "FGeometryInit.hh"
  37 #include "G4VPhysicalVolume.hh"
  38 #include "FluggNavigator.hh"
  39 #include "G4ThreeVector.hh"
  40 #include "globals.hh"
  41
  42
  43 void nrmlwr(G4double& /*pSx*/, G4double& /*pSy*/, G4double& /*pSz*/,
  44             G4double& pVx, G4double& pVy, G4double& pVz,
  45             G4double* norml, const G4int& oldReg,
  46             const G4int& /*newReg*/, G4int& flagErr)
  47 {
  48   //flag
  49 #ifdef G4GEOMETRY_DEBUG
  50   G4cout << "============ NRMLWR-DBG =============" << G4endl;
  51 #endif
  52
  53   //dummy variables
  54   flagErr=0;
  55
  56   //navigator pointer
  57   static FGeometryInit * ptrGeoInit;
  58   ptrGeoInit = FGeometryInit::GetInstance();
  59   FluggNavigator * ptrNavig = ptrGeoInit->getNavigatorForTracking();
  60
  61   //variables
  62   G4ThreeVector normalLoc;
  63   G4ThreeVector normalGlob;
  64
  65   //normal computing
  66   G4bool valid;
  67   normalLoc=ptrNavig->GetLocalExitNormal(&valid);
  68   if(valid) {
  69     normalLoc *= -1;
  70
  71     //global cooordinates normal
  72     normalGlob = ptrNavig->GetLocalToGlobalTransform().
  73       TransformAxis(normalLoc);
  74   }
  75   else {
  76     G4VPhysicalVolume *touchvolume;
  77     G4ThreeVector theLocalPoint;
  78
  79     if(ptrNavig->EnteredDaughterVolume()) {
  80       //volume from touchable history
  81       G4TouchableHistory *ptrTouchableHistory=ptrGeoInit->
  82         GetTouchableHistory();
  83       touchvolume = ptrTouchableHistory->GetVolume();
  84
  85       //local point from navigator and normal
  86       theLocalPoint = ptrNavig->GetCurrentLocalCoordinate();
  87       normalLoc = touchvolume->GetLogicalVolume()->GetSolid()->
  88         SurfaceNormal(theLocalPoint);
  89
  90       //global cooordinates normal
  91       normalGlob = ptrNavig->GetLocalToGlobalTransform().
  92         TransformAxis(normalLoc);
  93     }
  94     else {
  95       //volume from old history
  96       const G4NavigationHistory * ptrOldNavHist =
  97         ptrGeoInit->GetOldNavHist()->GetHistory();
  98       touchvolume = ptrOldNavHist->GetTopVolume();
  99
 100       if(!touchvolume) {
 101         // Old history has been reseted by LOOKZ relocation,
 102         // so is necessary to track back and forward to find
 103         // the right histories.
 104
 105
 106         ////////////  COMPUTE STEP BACKWARD //////////////////
 107         G4ThreeVector theGlobalPoint = ptrNavig->GetLocalToGlobalTransform().
 108           TransformPoint(ptrNavig->GetCurrentLocalCoordinate());
 109
 110         //compute step and new location
 111         G4ThreeVector pVec(pVx,pVy,pVz);
 112         pVec=-pVec;
 113         G4double appStep = 0;
 114         G4double safety = 0;
 115         G4bool onBoundary = false;
 116         G4double physStep = 1000000000;
 117         G4int newRegStep;
 118         G4ThreeVector partLoc =  theGlobalPoint;
 119         G4bool fErr=false;
 120
 121 #ifdef G4GEOMETRY_DEBUG
 122         G4cout << "Old history not found" << G4endl;
 123         G4cout << "* NRML needs boundary-crossing: computing step backward..."
 124                << G4endl;
 125 #endif
 126
 127         //compute step and location
 128         newRegStep=StepAndLocation(partLoc,pVec,physStep,
 129                                    appStep,safety,onBoundary,
 130                                    fErr,oldReg);
 131
 132         G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::
 133           GetInstance();
 134
 135         if(appStep<physStep && newRegStep!=G4int(pVolStore->size())+1) {
 136           //end-step point is on boundary between volumes;
 137           //==> update step-histories
 138
 139 #ifdef G4GEOMETRY_DEBUG
 140           G4cout << "* updating step-histories" << G4endl;
 141 #endif
 142
 143           G4TouchableHistory * ptrTouchHist =
 144             ptrGeoInit->GetTouchableHistory();
 145           ptrGeoInit->UpdateHistories(ptrTouchHist->GetHistory(),2);
 146         }
 147         else {
 148 #ifdef G4GEOMETRY_DEBUG
 149           G4cout << "ERROR! Boundary not found" << G4endl;
 150 #endif
 151         }
 152
 153 #ifdef G4GEOMETRY_DEBUG
 154         G4cout << "* computing step forward..." << G4endl;
 155 #endif
 156         pVec=-pVec;
 157         safety = 0;
 158         onBoundary = false;
 159
 160         //compute step and location for boundary crossing
 161         newRegStep=StepAndLocation(partLoc,pVec,physStep,
 162                                    appStep,safety,onBoundary,fErr,oldReg);
 163         if(appStep<physStep) {
 164           //end-step point is on boundary between volumes;
 165           //==> update step-histories
 166
 167 #ifdef G4GEOMETRY_DEBUG
 168           G4cout << "* updating step-histories" << G4endl;
 169 #endif
 170
 171           G4TouchableHistory * ptrTouchHist =
 172             ptrGeoInit->GetTouchableHistory();
 173           ptrGeoInit->UpdateHistories(ptrTouchHist->GetHistory(),2);
 174         }
 175       }
 176
 177       // now touchvolume exist.
 178       // local point from navigator and global point
 179       // N.B. if particle has exited world volume,
 180       // FluggNavigator doesn't update lastLocatedPoint.
 181       // So be carefull in building geometry always to have a
 182       // big world volume that fluka won't exit.
 183
 184       touchvolume = ptrOldNavHist->GetTopVolume();
 185
 186       G4ThreeVector theGlobalPoint = ptrNavig->GetLocalToGlobalTransform().
 187         TransformPoint(ptrNavig->GetCurrentLocalCoordinate());
 188       theLocalPoint = ptrOldNavHist->GetTopTransform().
 189         TransformPoint(theGlobalPoint);
 190       normalLoc = (touchvolume->GetLogicalVolume()->GetSolid()->
 191                    SurfaceNormal(theLocalPoint));
 192       normalLoc *= -1;
 193
 194       //global cooordinates normal
 195       normalGlob = ptrOldNavHist->GetTopTransform().
 196         Inverse().TransformAxis(normalLoc);
 197     }
 198   }
 199
 200   //return normal:
 201   norml[0]=G4double(normalGlob.x());
 202   norml[1]=G4double(normalGlob.y());
 203   norml[2]=G4double(normalGlob.z());
 204
 205   //for debugging
 206 #ifdef G4GEOMETRY_DEBUG
 207   G4cout << "Normal: " << normalGlob << G4endl;
 208 #endif
 209 }
 210
 211