]>
Commit | Line | Data |
---|---|---|
26911512 | 1 | #include "WrapUtils.hh" |
26911512 | 2 | #include "FGeometryInit.hh" |
26911512 | 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 | //////////////////////////////////////////////////////////////////////// | |
0edf14e7 | 172 | std::ostream& PrintHeader(std::ostream& os, const char* title) { |
26911512 | 173 | os << "*\n" << "*\n" << "*\n"; |
174 | os << "********************* " << title << " *********************\n" | |
175 | << "*\n"; | |
176 | os << "*...+....1....+....2....+....3....+....4....+....5....+....6....+....7..." | |
36c70081 | 177 | << G4endl; |
178 | os << "*" << G4endl; | |
26911512 | 179 | |
180 | return os; | |
181 | } | |
182 | ||
183 | //////////////////////////////////////////////////////////////////////// | |
26d97e06 | 184 | // ToFlukaString converts an string to something usefull in FLUKA: |
185 | // * Capital letters | |
186 | // * Only 8 letters | |
187 | // * Replace ' ' by '_' | |
26911512 | 188 | //////////////////////////////////////////////////////////////////////// |
26d97e06 | 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); | |
26911512 | 195 | } |
26d97e06 | 196 | return str; |
26911512 | 197 | } |