Flugg/WrapUtils.cxx

   1 #include "WrapUtils.hh"
   2 #include "FGeometryInit.hh"
   3
   4 ////////////////////////////////////////////////////////////////////////
   5 // StepAndLocation
   6 ////////////////////////////////////////////////////////////////////////
   7
   8 G4int StepAndLocation(G4ThreeVector &partLoc, const G4ThreeVector &pVec,
   9                       const G4double &physStep, G4double &appStep,
  10                       G4double &safety, G4bool & onBound,
  11                       G4bool & flagErr, const G4int & oldReg)
  12 {
  13   //NB onBound=true in the particular case when particle is in boundary
  14   //rounding BUT geometry hasn't crossed it.
  15   //flagErr=true when particle is on boundary but step=infinity,
  16   //newSafety>10.e-10 and lLocate function locates end step point in
  17   // a new region. This is absurb!
  18
  19
  20   //Geoinit, Navigator, TouchableHistory, etc. pointers
  21   static FGeometryInit * ptrGeoInit = FGeometryInit::GetInstance();
  22   FluggNavigator * ptrNavig = ptrGeoInit->getNavigatorForTracking();
  23   G4TouchableHistory * ptrTouchableHistory = ptrGeoInit->GetTouchableHistory();
  24   G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
  25
  26   //compute step
  27   G4double step = ptrNavig->ComputeStep(partLoc,pVec,physStep,safety);
  28
  29   //compute approved step
  30   if(step<physStep) {
  31     //NB:to check if point is really on boundary:
  32     //see G4AuxiliaryNavServices.hh
  33 #ifdef G4GEOMETRY_DEBUG
  34     G4cout << "* Boundary crossing!" << G4endl;
  35 #endif
  36
  37     appStep=step;
  38     ptrNavig->SetGeometricallyLimitedStep();
  39
  40     //if SetGeometricallyLimitedStep() is called, next
  41     //LocateGlobalPointAndUpdateTouchable(...,true) will
  42     //start searching from the end of the previous step,
  43     //otherwise will start from last located point.
  44     //(e.g.PropagatorInField and SteppingManger.
  45   }
  46   else {
  47     //step=kInfinity (the nearest boundary is >physStep)
  48     //Fluka always wants step lenght approved, so:
  49     appStep=physStep;
  50   }
  51
  52   //new position
  53   partLoc+=(pVec*appStep);
  54
  55   //locate point and update touchable history
  56   ptrNavig->LocateGlobalPointAndUpdateTouchable(partLoc,0,
  57                                                 ptrTouchableHistory,true);
  58
  59   G4VPhysicalVolume * touchvolume=ptrTouchableHistory->GetVolume();
  60
  61   //if volume not found, out of mother volume:
  62   //returns [number of volumes]+1
  63   G4int newReg;
  64   if(!touchvolume) {
  65 #ifdef G4GEOMETRY_DEBUG
  66     G4cout << "Out of mother volume! (G1WR)"  << G4endl;
  67 #endif
  68     G4int numVol = G4int(pVolStore->size());
  69     newReg = numVol + 1;
  70   }
  71   else {
  72 #ifdef G4GEOMETRY_DEBUG
  73     G4cout <<"Volume after step: " <<touchvolume->GetName() << G4endl;
  74 #endif
  75
  76     //compute new volume index
  77     G4int volIndex = ~0;
  78     for (unsigned int i=0; i<pVolStore->size(); i++)
  79       if ((*pVolStore)[i] == touchvolume)
  80         volIndex = i;
  81     if (volIndex==(~0)) {
  82       G4cerr << "ERROR in FLUGG: Problem in routine WrapG1 trying to find volume after step" << G4endl;
  83       exit(-999);
  84     }
  85     newReg=volIndex+1;
  86   }
  87
  88   onBound=false;
  89   flagErr=false;
  90   //check if position is in a boundary rounding while volume isn't changed
  91   if(step==kInfinity || step==physStep) {
  92     G4ThreeVector globalpoint = ptrNavig->GetLocalToGlobalTransform().
  93       TransformPoint(ptrNavig->GetCurrentLocalCoordinate());
  94
  95     //compute new safety
  96     G4double newSafetydbg = ptrNavig->ComputeSafety(globalpoint,DBL_MAX );
  97
  98 #ifdef G4GEOMETRY_DEBUG
  99     G4cout << "|From StepAndLocation:"  << G4endl;
 100     G4cout << "|step=" << step << G4endl;
 101     G4cout << "|phystep=" << physStep << G4endl;
 102     G4cout << "|safety=" << safety << G4endl;
 103     G4cout << "|newSafetydbg =" << newSafetydbg << G4endl;
 104 #endif
 105
 106     if(newSafetydbg<1.0e-10)
 107       onBound=true;
 108     else if(newReg != oldReg) {
 109 #ifdef G4GEOMETRY_DEBUG
 110       G4cout << "New Volume but ComputeStep didn't notice!" << G4endl;
 111 #endif
 112       flagErr=true;
 113     }
 114   }
 115
 116   //return
 117   return newReg;
 118 }
 119
 120
 121
 122
 123 ////////////////////////////////////////////////////////////////////////
 124 // EqualHistories
 125 ////////////////////////////////////////////////////////////////////////
 126
 127 bool EqualHistories(const G4NavigationHistory* ptrFirstHist,
 128                     const G4NavigationHistory* ptrSecHist)
 129 {
 130 #ifdef G4GEOMETRY_DEBUG
 131   G4cout << "* Testing if histories are equal..."  << G4endl;
 132 #endif
 133
 134   G4int depth1 = ptrFirstHist->GetDepth();
 135   G4int depth2 = ptrSecHist->GetDepth();
 136
 137   if(depth1!=depth2)
 138     return false;
 139
 140   for(G4int w=0;w<=depth1;w++) {
 141     if (*(ptrFirstHist->GetVolume(w))==*(ptrSecHist->GetVolume(w))) {
 142       //kNormal volume
 143       if(ptrFirstHist->GetVolumeType(w) == kNormal) {
 144         if ( ptrFirstHist->GetVolume(w)->GetCopyNo() !=
 145              ptrSecHist->GetVolume(w)->GetCopyNo() )
 146           return false;
 147       }
 148
 149       //Replica or Parametric volume
 150       else {
 151         if ( ptrFirstHist->GetReplicaNo(w) !=
 152              ptrSecHist->GetReplicaNo(w) )
 153           return false;
 154       }
 155     }
 156     else
 157       return false;
 158   }
 159
 160 #ifdef G4GEOMETRY_DEBUG
 161   G4cout << "Histories are equal!" << G4endl;
 162 #endif
 163
 164   return true;
 165 }
 166
 167
 168
 169 ////////////////////////////////////////////////////////////////////////
 170 // PrintHeader
 171 ////////////////////////////////////////////////////////////////////////
 172 std::ostream& PrintHeader(std::ostream& os, const char* title) {
 173   os << "*\n" << "*\n" << "*\n";
 174   os << "*********************  " << title << " *********************\n"
 175      << "*\n";
 176   os << "*...+....1....+....2....+....3....+....4....+....5....+....6....+....7..."
 177      << G4endl;
 178   os << "*" << G4endl;
 179
 180   return os;
 181 }
 182
 183 ////////////////////////////////////////////////////////////////////////
 184 // ToFlukaString converts an string to something usefull in FLUKA:
 185 // * Capital letters
 186 // * Only 8 letters
 187 // * Replace ' ' by '_'
 188 ////////////////////////////////////////////////////////////////////////
 189 G4String ToFlukaString(G4String str) {
 190   str = str.substr(0,8);
 191   str.toUpper();
 192   for (unsigned int i = 0; i < str.length(); i++) {
 193     if (str.at(i) == ' ')
 194       str.replace(i,1,"_",1);
 195   }
 196   return str;
 197 }