1 #include "WrapUtils.hh"
2 #include "FGeometryInit.hh"
4 ////////////////////////////////////////////////////////////////////////
6 ////////////////////////////////////////////////////////////////////////
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)
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!
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();
27 G4double step = ptrNavig->ComputeStep(partLoc,pVec,physStep,safety);
29 //compute approved step
31 //NB:to check if point is really on boundary:
32 //see G4AuxiliaryNavServices.hh
33 #ifdef G4GEOMETRY_DEBUG
34 G4cout << "* Boundary crossing!" << G4endl;
38 ptrNavig->SetGeometricallyLimitedStep();
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.
47 //step=kInfinity (the nearest boundary is >physStep)
48 //Fluka always wants step lenght approved, so:
53 partLoc+=(pVec*appStep);
55 //locate point and update touchable history
56 ptrNavig->LocateGlobalPointAndUpdateTouchable(partLoc,0,
57 ptrTouchableHistory,true);
59 G4VPhysicalVolume * touchvolume=ptrTouchableHistory->GetVolume();
61 //if volume not found, out of mother volume:
62 //returns [number of volumes]+1
65 #ifdef G4GEOMETRY_DEBUG
66 G4cout << "Out of mother volume! (G1WR)" << G4endl;
68 G4int numVol = G4int(pVolStore->size());
72 #ifdef G4GEOMETRY_DEBUG
73 G4cout <<"Volume after step: " <<touchvolume->GetName() << G4endl;
76 //compute new volume index
78 for (unsigned int i=0; i<pVolStore->size(); i++)
79 if ((*pVolStore)[i] == touchvolume)
82 G4cerr << "ERROR in FLUGG: Problem in routine WrapG1 trying to find volume after step" << G4endl;
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());
96 G4double newSafetydbg = ptrNavig->ComputeSafety(globalpoint,DBL_MAX );
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;
106 if(newSafetydbg<1.0e-10)
108 else if(newReg != oldReg) {
109 #ifdef G4GEOMETRY_DEBUG
110 G4cout << "New Volume but ComputeStep didn't notice!" << G4endl;
123 ////////////////////////////////////////////////////////////////////////
125 ////////////////////////////////////////////////////////////////////////
127 bool EqualHistories(const G4NavigationHistory* ptrFirstHist,
128 const G4NavigationHistory* ptrSecHist)
130 #ifdef G4GEOMETRY_DEBUG
131 G4cout << "* Testing if histories are equal..." << G4endl;
134 G4int depth1 = ptrFirstHist->GetDepth();
135 G4int depth2 = ptrSecHist->GetDepth();
140 for(G4int w=0;w<=depth1;w++) {
141 if (*(ptrFirstHist->GetVolume(w))==*(ptrSecHist->GetVolume(w))) {
143 if(ptrFirstHist->GetVolumeType(w) == kNormal) {
144 if ( ptrFirstHist->GetVolume(w)->GetCopyNo() !=
145 ptrSecHist->GetVolume(w)->GetCopyNo() )
149 //Replica or Parametric volume
151 if ( ptrFirstHist->GetReplicaNo(w) !=
152 ptrSecHist->GetReplicaNo(w) )
160 #ifdef G4GEOMETRY_DEBUG
161 G4cout << "Histories are equal!" << G4endl;
169 ////////////////////////////////////////////////////////////////////////
171 ////////////////////////////////////////////////////////////////////////
172 std::ostream& PrintHeader(std::ostream& os, const char* title) {
173 os << "*\n" << "*\n" << "*\n";
174 os << "********************* " << title << " *********************\n"
176 os << "*...+....1....+....2....+....3....+....4....+....5....+....6....+....7..."
183 ////////////////////////////////////////////////////////////////////////
184 // ToFlukaString converts an string to something usefull in FLUKA:
187 // * Replace ' ' by '_'
188 ////////////////////////////////////////////////////////////////////////
189 G4String ToFlukaString(G4String str) {
190 str = str.substr(0,8);
192 for (unsigned int i = 0; i < str.length(); i++) {
193 if (str.at(i) == ' ')
194 str.replace(i,1,"_",1);