Flugg/WrapLookMG.cxx

   1
   2 // Flugg tag
   3
   4 ///////////////////////////////////////////////////////////////////
   5 //
   6 // WrapLookMG.hh - Sara Vanini 26/X/99
   7 //
   8 // Wrapper for localisation of particle for magnetic field tracking
   9 //
  10 // modified 13/IV/00: check if point belongs to one of the lattice
  11 // histories stored in jrLtGeant
  12 //
  13 // modified 24.10.01: by I. Hrivnacova
  14 //   functions declarations separated from implementation
  15 //   (moved to Wrappers.hh);
  16 // modified 17/06/02: by I. Gonzalez. STL migration
  17 //
  18 ///////////////////////////////////////////////////////////////////
  19
  20
  21 #include "Wrappers.hh"
  22 #include "FGeometryInit.hh"
  23 #include "NavHistWithCount.hh"
  24 #include "G4PhysicalVolumeStore.hh"
  25 #include "G4NormalNavigation.hh"
  26 #include "G4VoxelNavigation.hh"
  27 #include "G4ParameterisedNavigation.hh"
  28 #include "G4ReplicaNavigation.hh"
  29 #include "globals.hh"
  30
  31
  32 //auxiliary function declarations
  33 G4bool PointLocate(const G4NavigationHistory &,
  34                           const G4ThreeVector &,
  35                           const G4ThreeVector *,
  36                           G4int &);
  37 EVolume CharacteriseDaughters(const G4LogicalVolume *);
  38
  39
  40 void lkmgwr(G4double& pSx, G4double& pSy, G4double& pSz,
  41             G4double* pV, const G4int& oldReg, const G4int& oldLttc,
  42             G4int& flagErr, G4int& newReg, G4int& newLttc)
  43 {
  44 //flag
  45 #ifdef G4GEOMETRY_DEBUG
  46     G4cout<<"============= LKMGWR =============="<<G4endl;
  47 #endif
  48
  49     //Geoinit, Navigator, etc. pointers
  50     static FGeometryInit * ptrGeoInit = FGeometryInit::GetInstance();
  51     G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
  52
  53     //setting variables (and dimension: Fluka uses cm.!)
  54     G4ThreeVector globalPointcm(pSx,pSy,pSz);
  55     G4ThreeVector globalPoint =  globalPointcm * 10.; //in mm.
  56     G4ThreeVector globalDirection(pV[0],pV[1],pV[2]);
  57
  58     //get jrLtGeant array and initialize variables
  59     G4int * jrLtGeant = ptrGeoInit->GetJrLtGeantArray();
  60     G4int LttcFlagGeant = ptrGeoInit->GetLttcFlagGeant();
  61     G4bool belongToVolume = false;
  62     G4int i = LttcFlagGeant;
  63     G4int regionIn = 0;
  64     newReg = -1;
  65     newLttc = -1;
  66
  67     while(!belongToVolume && i>=0) {
  68       //get history from jrLtGeant
  69       NavHistWithCount * ptrNavHistCount =
  70         reinterpret_cast<NavHistWithCount*>(jrLtGeant[i]);
  71       const G4NavigationHistory * ptrNavHist =
  72         ptrNavHistCount->GetNavHistPtr();
  73
  74       //check if globalPoint belongs to volume (can't call
  75       //LocateGlobalPoint... because of flag settings)
  76       belongToVolume = PointLocate(*ptrNavHist,globalPoint,
  77                                    &globalDirection,regionIn);
  78
  79
  80       //if point belongs to surface, check scalar product direction*normal
  81       if(regionIn==-100) {
  82 #ifdef G4GEOMETRY_DEBUG
  83         G4cout<<"On surface!"<<G4endl;
  84 #endif
  85
  86         //if entering, then point belongs to volume,
  87         //if exiting, point doesn't belong to volume.
  88         G4int oldReg,flagErr,reg;
  89         G4double x,y,z,px,py,pz;
  90         G4double * norml = new G4double[3];
  91         x = globalPoint.x();
  92         y = globalPoint.y();
  93         z = globalPoint.z();
  94         px = globalDirection.x();
  95         py = globalDirection.y();
  96         pz = globalDirection.z();
  97
  98         nrmlwr(x,y,z,px,py,pz,norml,oldReg,reg,flagErr);
  99
 100         G4ThreeVector normal(norml[0],norml[1],norml[2]);
 101 #ifdef G4GEOMETRY_DEBUG
 102         // G4cout<<"Scalar product="<<globalDirection.dot(normal)<<G4endl;
 103 #endif
 104         if(globalDirection.dot(normal)>0) {
 105 #ifdef G4GEOMETRY_DEBUG
 106           G4cout<<"entering volume!"<<G4endl;
 107 #endif
 108           G4int volIndex = ~0;
 109           for (unsigned int i=0; i<pVolStore->size(); i++)
 110             if ((*pVolStore)[i] == ptrNavHist->GetTopVolume())
 111               volIndex = i;
 112           if (volIndex==(~0)) {
 113             G4cerr << "FLUGG: Problem in routine WrapG1 tryingto find volume after step" << G4endl;
 114             exit(-999);
 115           }
 116           //regionIn = G4int(pVolStore->index(ptrNavHist->GetTopVolume()))+1;
 117           regionIn = volIndex+1;
 118           belongToVolume = true;
 119         }
 120       }
 121
 122       i -= 1;
 123     }
 124
 125     //output variables
 126     if(belongToVolume) {
 127       newReg = regionIn;
 128       newLttc = jrLtGeant[i+1];
 129     }
 130
 131     flagErr = pVolStore->size() + 2;      //flagErr=fluka region number + 1
 132
 133
 134 #ifdef G4GEOMETRY_DEBUG
 135     G4cout << "Global point (cm): " << globalPointcm << G4endl;
 136     G4cout << "Direction: " << globalDirection << G4endl;
 137     if(newReg!=-1)
 138       G4cout << "Point belongs to region " << newReg
 139              << " ptr history=" << newLttc << G4endl;
 140     else
 141       G4cout << "No containing volume found!" << G4endl;
 142 #endif
 143 }
 144
 145
 146
 147
 148 //returns false if the point doesn't belong to history top volume (can either
 149 //belong to one of its daughter volumes or none), otherwise returns true.
 150
 151 G4bool PointLocate(const G4NavigationHistory & blockedHistConst,
 152                    const G4ThreeVector & globalPoint,
 153                    const G4ThreeVector * pGlobalDirection,
 154                    G4int & reg)
 155 {
 156
 157   //variables
 158   G4NavigationHistory * fHistory = new G4NavigationHistory(blockedHistConst);
 159   G4bool belongVolume;
 160
 161   //G4 flags resetted (see: ResetStackAndState)
 162   G4bool fExiting = false;
 163   G4VPhysicalVolume * fBlockedPhysicalVolume=0;
 164   G4int fBlockedReplicaNo=-1;
 165   G4bool fLocatedOnEdge=false;
 166
 167   G4bool notKnownContained=true,noResult;
 168   G4VPhysicalVolume *targetPhysical;
 169   G4LogicalVolume *targetLogical;
 170   G4VSolid *targetSolid;
 171   G4ThreeVector localPoint,localDirection;
 172   EInside insideCode;
 173
 174
 175   //Helpers/Utility classes
 176   G4NormalNavigation  fnormalNav;
 177   G4VoxelNavigation fvoxelNav;
 178   G4ParameterisedNavigation fparamNav;
 179   G4ReplicaNavigation freplicaNav;
 180   G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
 181
 182   //local variables
 183   localPoint=fHistory->GetTopTransform().TransformPoint(globalPoint);
 184   localDirection=fHistory->GetTopTransform().TransformPoint(*pGlobalDirection);
 185
 186   //search if volume contains point
 187   if (fHistory->GetTopVolumeType()!=kReplica) {
 188     targetSolid=fHistory->GetTopVolume()->GetLogicalVolume()->GetSolid();
 189     insideCode=targetSolid->Inside(localPoint);
 190   }
 191   else {
 192     insideCode=freplicaNav.BackLocate(*fHistory,globalPoint,
 193                                       localPoint,fExiting,notKnownContained);
 194     // !CARE! if notKnownContained returns false then the point is within
 195     // the containing placement volume of the replica(s). If insidecode
 196     // will result in the history being backed up one level, then the
 197     // local point returned is the point in the system of this new level
 198 #ifdef G4GEOMETRY_DEBUG
 199     G4cout << "Replica: fExiting=" << fExiting << G4endl;
 200     G4cout << "         notKnownContained=" << notKnownContained << G4endl;
 201 #endif
 202   }
 203
 204   if (insideCode==kOutside)
 205     belongVolume = false;
 206
 207   else if (insideCode==kSurface) {
 208     belongVolume = false;
 209     reg = -100;
 210   }
 211   else {
 212     // Search downwards in daughter volumes
 213     // until deepest containing volume found
 214     //
 215     // 3 Cases:
 216     //
 217     // o Parameterised daughters
 218     //   =>Must be one G4PVParameterised daughter & voxels
 219     // o Positioned daughters & voxels
 220     // o Positioned daughters & no voxels
 221
 222     belongVolume = true;
 223     noResult=true;
 224
 225     do {
 226       // Determine `type' of current mother volume
 227       targetPhysical=fHistory->GetTopVolume();
 228       targetLogical=targetPhysical->GetLogicalVolume();
 229       switch(CharacteriseDaughters(targetLogical)) {
 230       case kNormal:
 231         if (targetLogical->GetVoxelHeader()) {
 232           noResult=fvoxelNav.LevelLocate(*fHistory,
 233                                          fBlockedPhysicalVolume,
 234                                          fBlockedReplicaNo,
 235                                          globalPoint,
 236                                          pGlobalDirection,
 237                                          fLocatedOnEdge,
 238                                          localPoint);
 239         }
 240         else {
 241           noResult=fnormalNav.LevelLocate(*fHistory,
 242                                           fBlockedPhysicalVolume,
 243                                           fBlockedReplicaNo,
 244                                           globalPoint,
 245                                           pGlobalDirection,
 246                                           fLocatedOnEdge,
 247                                           localPoint);
 248         }
 249         break;
 250       case kReplica:
 251         noResult=freplicaNav.LevelLocate(*fHistory,
 252                                          fBlockedPhysicalVolume,
 253                                          fBlockedReplicaNo,
 254                                          globalPoint,
 255                                          pGlobalDirection,
 256                                          fLocatedOnEdge,
 257                                          localPoint);
 258         break;
 259       case kParameterised:
 260         noResult=fparamNav.LevelLocate(*fHistory,
 261                                        fBlockedPhysicalVolume,
 262                                        fBlockedReplicaNo,
 263                                        globalPoint,
 264                                        pGlobalDirection,
 265                                        fLocatedOnEdge,
 266                                        localPoint);
 267         break;
 268       } //switch
 269
 270       // LevelLocate search in the first daughter level.
 271       // LevelLocate returns noResult=true if it finds a daughter volume
 272       // in which globalPoint is inside (or on the surface). So point
 273       // doesn't belong only to mother volume ==> belongVolume=false.
 274
 275       if (noResult) {
 276         // The blocked volume no longer valid - it was for another level
 277         fBlockedPhysicalVolume= 0;
 278         fBlockedReplicaNo= -1;
 279         belongVolume=false;
 280       }
 281     } while (noResult);
 282
 283     //on exit targetPhysical is the volume globalPoint belongs to;
 284     G4int volIndex = ~0;
 285     for (unsigned int i=0; i<pVolStore->size(); i++)
 286       if ((*pVolStore)[i] == targetPhysical)
 287         volIndex = i;
 288     if (volIndex==(~0)) {
 289       G4cerr << "FLUGG: Problem in routine WrapG1 tryingto find volume after step" << G4endl;
 290       exit(-999);
 291     }
 292     //reg = G4int(pVolStore->index(targetPhysical))+1;
 293     reg = volIndex+1;
 294
 295   }
 296
 297   delete fHistory;
 298   return belongVolume;
 299 }
 300
 301 EVolume CharacteriseDaughters(const G4LogicalVolume *pLog)
 302 {
 303   EVolume type;
 304   EAxis axis;
 305   G4int nReplicas;
 306   G4double width,offset;
 307   G4bool consuming;
 308   G4VPhysicalVolume *pVol;
 309
 310   if (pLog->GetNoDaughters()==1) {
 311     pVol=pLog->GetDaughter(0);
 312     if (pVol->IsReplicated()) {
 313       pVol->GetReplicationData(axis,nReplicas,width,offset,consuming);
 314       type=(consuming) ? kReplica : kParameterised;
 315     }
 316     else
 317       type=kNormal;
 318   }
 319   else
 320     type=kNormal;
 321   return type;
 322 }
 323
 324
 325