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   if(ptrNavig->IsExitNormalValid()) {
  67     normalLoc=ptrNavig->GetLocalExitNormal();
  68     normalLoc *= -1;
  69
  70     //global cooordinates normal
  71     normalGlob = ptrNavig->GetLocalToGlobalTransform().
  72       TransformAxis(normalLoc);
  73   }
  74   else {
  75     G4VPhysicalVolume *touchvolume;
  76     G4ThreeVector theLocalPoint;
  77
  78     if(ptrNavig->EnteredDaughterVolume()) {
  79       //volume from touchable history
  80       G4TouchableHistory *ptrTouchableHistory=ptrGeoInit->
  81         GetTouchableHistory();
  82       touchvolume = ptrTouchableHistory->GetVolume();
  83
  84       //local point from navigator and normal
  85       theLocalPoint = ptrNavig->GetCurrentLocalCoordinate();
  86       normalLoc = touchvolume->GetLogicalVolume()->GetSolid()->
  87         SurfaceNormal(theLocalPoint);
  88
  89       //global cooordinates normal
  90       normalGlob = ptrNavig->GetLocalToGlobalTransform().
  91         TransformAxis(normalLoc);
  92     }
  93     else {
  94       //volume from old history
  95       const G4NavigationHistory * ptrOldNavHist =
  96         ptrGeoInit->GetOldNavHist()->GetHistory();
  97       touchvolume = ptrOldNavHist->GetTopVolume();
  98
  99       if(!touchvolume) {
 100         // Old history has been reseted by LOOKZ relocation,
 101         // so is necessary to track back and forward to find
 102         // the right histories.
 103
 104
 105         ////////////  COMPUTE STEP BACKWARD //////////////////
 106         G4ThreeVector theGlobalPoint = ptrNavig->GetLocalToGlobalTransform().
 107           TransformPoint(ptrNavig->GetCurrentLocalCoordinate());
 108
 109         //compute step and new location
 110         G4ThreeVector pVec(pVx,pVy,pVz);
 111         pVec=-pVec;
 112         G4double appStep = 0;
 113         G4double safety = 0;
 114         G4bool onBoundary = false;
 115         G4double physStep = 1000000000;
 116         G4int newRegStep;
 117         G4ThreeVector partLoc =  theGlobalPoint;
 118         G4bool fErr=false;
 119
 120 #ifdef G4GEOMETRY_DEBUG
 121         G4cout << "Old history not found" << G4endl;
 122         G4cout << "* NRML needs boundary-crossing: computing step backward..."
 123                << G4endl;
 124 #endif
 125
 126         //compute step and location
 127         newRegStep=StepAndLocation(partLoc,pVec,physStep,
 128                                    appStep,safety,onBoundary,
 129                                    fErr,oldReg);
 130
 131         G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::
 132           GetInstance();
 133
 134         if(appStep<physStep && newRegStep!=G4int(pVolStore->size())+1) {
 135           //end-step point is on boundary between volumes;
 136           //==> update step-histories
 137
 138 #ifdef G4GEOMETRY_DEBUG
 139           G4cout << "* updating step-histories" << G4endl;
 140 #endif
 141
 142           G4TouchableHistory * ptrTouchHist =
 143             ptrGeoInit->GetTouchableHistory();
 144           ptrGeoInit->UpdateHistories(ptrTouchHist->GetHistory(),2);
 145         }
 146         else {
 147 #ifdef G4GEOMETRY_DEBUG
 148           G4cout << "ERROR! Boundary not found" << G4endl;
 149 #endif
 150         }
 151
 152 #ifdef G4GEOMETRY_DEBUG
 153         G4cout << "* computing step forward..." << G4endl;
 154 #endif
 155         pVec=-pVec;
 156         safety = 0;
 157         onBoundary = false;
 158
 159         //compute step and location for boundary crossing
 160         newRegStep=StepAndLocation(partLoc,pVec,physStep,
 161                                    appStep,safety,onBoundary,fErr,oldReg);
 162         if(appStep<physStep) {
 163           //end-step point is on boundary between volumes;
 164           //==> update step-histories
 165
 166 #ifdef G4GEOMETRY_DEBUG
 167           G4cout << "* updating step-histories" << G4endl;
 168 #endif
 169
 170           G4TouchableHistory * ptrTouchHist =
 171             ptrGeoInit->GetTouchableHistory();
 172           ptrGeoInit->UpdateHistories(ptrTouchHist->GetHistory(),2);
 173         }
 174       }
 175
 176       // now touchvolume exist.
 177       // local point from navigator and global point
 178       // N.B. if particle has exited world volume,
 179       // FluggNavigator doesn't update lastLocatedPoint.
 180       // So be carefull in building geometry always to have a
 181       // big world volume that fluka won't exit.
 182
 183       touchvolume = ptrOldNavHist->GetTopVolume();
 184
 185       G4ThreeVector theGlobalPoint = ptrNavig->GetLocalToGlobalTransform().
 186         TransformPoint(ptrNavig->GetCurrentLocalCoordinate());
 187       theLocalPoint = ptrOldNavHist->GetTopTransform().
 188         TransformPoint(theGlobalPoint);
 189       normalLoc = (touchvolume->GetLogicalVolume()->GetSolid()->
 190                    SurfaceNormal(theLocalPoint));
 191       normalLoc *= -1;
 192
 193       //global cooordinates normal
 194       normalGlob = ptrOldNavHist->GetTopTransform().
 195         Inverse().TransformAxis(normalLoc);
 196     }
 197   }
 198
 199   //return normal:
 200   norml[0]=G4double(normalGlob.x());
 201   norml[1]=G4double(normalGlob.y());
 202   norml[2]=G4double(normalGlob.z());
 203
 204   //for debugging
 205 #ifdef G4GEOMETRY_DEBUG
 206   G4cout << "Normal: " << normalGlob << G4endl;
 207 #endif
 208 }
 209
 210