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 G4std::ostream& PrintHeader(G4std::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 // PrintMaterial
 185 ////////////////////////////////////////////////////////////////////////
 186 G4std::ostream& PrintMaterial(G4std::ostream& os, const char* title,
 187                        G4double Z, G4double A,
 188                        G4double density,
 189                        G4double index,
 190                        G4double N,
 191                        const char* name) {
 192
 193   os << setw10 << title;
 194
 195   os.setf(static_cast<G4std::ios::fmtflags>(0),G4std::ios::floatfield);
 196   if (Z < 0)
 197     os << setw10 << " ";
 198   else
 199     os << setw10
 200        << setfixed
 201        << G4std::setprecision(1)
 202        << Z;
 203
 204   if (A < 0)
 205     os << setw10 << " ";
 206   else
 207     os << setw10 << G4std::setprecision(3)
 208        << A;
 209
 210   os.setf(static_cast<G4std::ios::fmtflags>(0),G4std::ios::floatfield);
 211   os << setw10
 212      << setscientific
 213      << G4std::setprecision(3)
 214      << density;
 215
 216   os.setf(static_cast<G4std::ios::fmtflags>(0),G4std::ios::floatfield);
 217   os << setw10
 218      << setfixed
 219      << G4std::setprecision(1)
 220      << index;
 221
 222   os << setw10 << " ";
 223
 224   if (N < 0)
 225     os << setw10 << " ";
 226   else
 227     os << setw10 << N;
 228
 229   os << name << G4endl;
 230
 231   return os;
 232 }
 233
 234
 235 ////////////////////////////////////////////////////////////////////////
 236 // PrintCompund
 237 ////////////////////////////////////////////////////////////////////////
 238 G4std::ostream& PrintCompound(G4std::ostream& os, const char* title,
 239                        G4int count,
 240                        const char* name,
 241                        G4double fraction,
 242                        G4double index) {
 243
 244
 245   if(count==3) {
 246     os << name << G4endl;
 247     os << setw10 << "COMPOUND  ";
 248   }
 249
 250   os.setf(static_cast<G4std::ios::fmtflags>(0),G4std::ios::floatfield);
 251   os << setw10
 252      << setscientific
 253      << G4std::setprecision(2)
 254      << fraction;
 255   os.setf(static_cast<G4std::ios::fmtflags>(0),G4std::ios::floatfield);
 256   os << setw10
 257      << setfixed
 258      << G4std::setprecision(1)
 259      << index;
 260
 261   return os;
 262 }
 263