]>
Commit | Line | Data |
---|---|---|
26911512 | 1 | |
2 | // Flugg tag | |
3 | ||
4 | /////////////////////////////////////////////////////////////////// | |
5 | // | |
6 | // WrapG1.hh - Sara Vanini | |
7 | // | |
8 | // Wrapper for geometry tracking: returns approved step of | |
9 | // particle and alla variables that fluka G1 computes. | |
10 | // | |
11 | // modified 11/III/99 : included jrLtGeant array for storing | |
12 | // lattice histories; fixed jump-on-boundaries and back-scattering | |
13 | // modified 20/IV/99 : history counters of jrLt array are | |
14 | // incremented when created and decremented when check is required: | |
15 | // this for not deleting histories twice. | |
16 | // modified 18/X/99 : LocateGlobalPointWithinVolume used when position | |
17 | // is changed from last located point. | |
18 | // modified 1/III/00 : update utilities histories when crossing | |
19 | // identical volume boundaries. | |
20 | // | |
21 | // modified 22/III/00 : fixed LttcFlag and jrLt return values. | |
22 | // modified 12/VI/00 : end-step on Boundary bug fixed. | |
23 | // modified 5/VII/00 : boundary not seen by G4 geometry bug fixed. | |
24 | // modified 24.10.01: by I. Hrivnacova | |
25 | // functions declarations separated from implementation | |
26 | // (moved to Wrappers.hh); | |
27 | // | |
28 | ///////////////////////////////////////////////////////////////////// | |
29 | ||
30 | #include "Wrappers.hh" | |
31 | #include "WrapUtils.hh" | |
32 | #include "FGeometryInit.hh" | |
33 | #include "NavHistWithCount.hh" | |
34 | #include "G4VPhysicalVolume.hh" | |
35 | #include "G4NavigationHistory.hh" | |
36 | #include "FluggNavigator.hh" | |
37 | #include "G4PhysicalVolumeStore.hh" | |
38 | #include "globals.hh" | |
39 | ||
40 | ||
41 | void g1wr(G4double& pSx, G4double& pSy, G4double& pSz, G4double* pV, | |
42 | G4int& oldReg, const G4int& oldLttc, G4double& propStep, | |
43 | G4int& nascFlag, G4double& retStep, G4int& newReg, | |
44 | G4double& saf, G4int& newLttc, G4int& LttcFlag, | |
45 | G4double* sLt, G4int* jrLt) | |
46 | { | |
47 | //flag | |
48 | #ifdef G4GEOMETRY_DEBUG | |
49 | G4cout<<"============= G1WR =============="<<G4endl; | |
50 | #endif | |
51 | ||
52 | //static G4int count=0; | |
53 | //count+=1; | |
54 | //G4cout<<"contatore G1="<<count<<G4endl; | |
55 | ||
56 | ///////////////////////// INPUT /////////////////////////// | |
57 | //Geoinit, Navigator, TouchableHistory pointers | |
58 | static FGeometryInit * ptrGeoInit = FGeometryInit::GetInstance(); | |
59 | FluggNavigator * ptrNavig = ptrGeoInit->getNavigatorForTracking(); | |
60 | G4TouchableHistory * ptrTouchHist = ptrGeoInit->GetTouchableHistory(); | |
61 | G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance(); | |
62 | ||
63 | //setting variables (and dimension: Fluka uses cm.!) | |
64 | G4ThreeVector partLoc(pSx,pSy,pSz); | |
65 | partLoc *= 10.0; // in millimeters! | |
66 | static G4ThreeVector partLocOld = partLoc; | |
67 | static G4ThreeVector oldLocalPoint = | |
68 | ptrNavig->ComputeLocalPoint(partLocOld); | |
69 | ||
70 | G4ThreeVector pVec(pV[0],pV[1],pV[2]); | |
71 | const G4double physStep=G4double(propStep*10.); | |
72 | ||
73 | G4int LttcFlagGeant = ptrGeoInit->GetLttcFlagGeant(); | |
74 | ||
75 | ||
76 | #ifdef G4GEOMETRY_DEBUG | |
77 | G4cout.precision(10); | |
78 | G4cout << "Position (cm):" << pSx << "," << pSy << "," << pSz << G4endl; | |
79 | G4cout << "Direction: " << pVec << G4endl; | |
80 | G4cout << "Proposed step :"<< propStep << G4endl; | |
81 | #endif | |
82 | ||
83 | ||
84 | ||
85 | ///////////////////// FLUKA/G4 REGIONS COMPARISON ////////////////////// | |
86 | //get oldReg pointer and G4volume of previous step end-point pointer | |
87 | G4VPhysicalVolume * ptrOldReg = (*pVolStore)[oldReg-1]; | |
88 | G4VPhysicalVolume * lastVolume = ptrTouchHist->GetVolume(); | |
89 | G4int touVolCpNb = 0; | |
90 | if (lastVolume) | |
91 | touVolCpNb = lastVolume->GetCopyNo(); | |
92 | ||
93 | #ifdef G4GEOMETRY_DEBUG | |
94 | G4int fluVolCpNb = ptrOldReg->GetCopyNo(); | |
95 | G4cout << "Fluka volume before step: " << ptrOldReg->GetName() | |
96 | << "," << fluVolCpNb<< G4endl; | |
97 | if(lastVolume) | |
98 | G4cout << "G4 Touch Hist volume: " | |
99 | << lastVolume->GetName() << "," << touVolCpNb << G4endl; | |
100 | G4cout << "------------------------------------------------" << G4endl; | |
101 | #endif | |
102 | ||
103 | ||
104 | //if volume is changed, this is a new particle tracking, or fluka tries | |
105 | //to reach a boundary more softly with the same particle. In this case | |
106 | //fluka restart tracking from old history, in general. For tracking in | |
107 | //lattice-volume fluka goes back to one of the previous lattice volumes. | |
108 | //Check if ptrOldReg is equal to old history, or to one of the N lattice | |
109 | //volumes stored in jrLt. Then reinitialise step-histories! Otherwise | |
110 | //must relocate. | |
111 | //NB. jrLtGeant stores lattice volume histories until LttcFlag==-1, | |
112 | //then all histories are checked and deleted, and array is reinitialised | |
113 | ||
114 | ||
115 | G4int haveHistNb = -1; | |
116 | G4int newRegErr=0; | |
117 | G4int indHist = LttcFlagGeant; | |
118 | G4int * jrLtGeant = ptrGeoInit->GetJrLtGeantArray(); | |
119 | ||
120 | ||
121 | G4NavigationHistory* ptrLttcHist=0; | |
122 | if(oldLttc) | |
123 | ptrLttcHist = reinterpret_cast<NavHistWithCount*>(oldLttc)->GetNavHistPtr(); | |
124 | ||
125 | ||
126 | while(indHist>=0 && haveHistNb==-1) { | |
127 | #ifdef G4GEOMETRY_DEBUG | |
128 | G4cout<<"Searching in jrLtc...."<<G4endl; | |
129 | #endif | |
130 | if(oldLttc==jrLtGeant[indHist]) | |
131 | haveHistNb=indHist; | |
132 | indHist-=1; | |
133 | } | |
134 | ||
135 | if(haveHistNb!=-1) { | |
136 | //fluka history found in jrLtGeant | |
137 | if(haveHistNb<LttcFlagGeant) { | |
138 | #ifdef G4GEOMETRY_DEBUG | |
139 | G4cout<<"* Fluka reaches boundary more softly..."<<G4endl; | |
140 | G4cout<<"* Re-initializing FluggNavigator history"<<G4endl; | |
141 | G4cout<<"* and updating step-histories"<<G4endl; | |
142 | #endif | |
143 | ||
144 | ptrNavig->UpdateNavigatorHistory(ptrLttcHist); | |
145 | ptrGeoInit->UpdateHistories(ptrLttcHist,0); | |
146 | } | |
147 | #ifdef G4GEOMETRY_DEBUG | |
148 | if(haveHistNb==LttcFlagGeant) | |
149 | G4cout << "Continuing step...." << G4endl; | |
150 | #endif | |
151 | jrLt[0]=oldLttc; | |
152 | } | |
153 | else { | |
154 | //not found fluka history in jrLttGeant! | |
155 | G4cout << "* ERROR! Geometry not correctly initialised in fluka history!" | |
156 | << G4endl; | |
157 | ||
158 | //relocation! | |
159 | ptrNavig->LocateGlobalPointAndUpdateTouchable(partLoc,0,ptrTouchHist,true); | |
160 | ||
161 | G4cout << "* ATTENTION: point relocation in: " | |
162 | << ptrTouchHist->GetVolume()->GetName() << G4endl; | |
163 | ||
164 | ptrGeoInit->UpdateHistories(ptrTouchHist->GetHistory(),1); | |
165 | ||
166 | if(ptrTouchHist->GetVolume() != ptrOldReg) { | |
167 | G4cout << "* ERROR! Point not in fluka volume!" << G4endl; | |
168 | newRegErr=-3; | |
169 | } | |
170 | ||
171 | //save new history in jrLt[0] and increment its counter | |
172 | #ifdef G4GEOMETRY_DEBUG | |
173 | G4cout << "* ISVHWR call to store new NavHistWithCount in jrLt[0]" | |
174 | << G4endl; | |
175 | #endif | |
176 | jrLt[0]=isvhwr(0,0); | |
177 | ||
178 | #ifdef G4GEOMETRY_DEBUG | |
179 | G4cout << "* CONHWR call to increment counter" << G4endl; | |
180 | #endif | |
181 | G4int incrCount2=1; | |
182 | conhwr(jrLt[0],&incrCount2); | |
183 | } | |
184 | ||
185 | ||
186 | //jrLtGeant - history check: decrement counter and delete | |
187 | //histories, if LttcFlag=-1, then reinitialise array with -1. | |
188 | if(LttcFlag==-1 && LttcFlagGeant>=0) { | |
189 | #ifdef G4GEOMETRY_DEBUG | |
190 | G4cout << "* CONHWR call to check and delete histories in jrLtGeant[]:" | |
191 | << G4endl; | |
192 | #endif | |
193 | ||
194 | for(G4int ind=0; ind<=LttcFlagGeant; ind++) { | |
195 | G4int incrCount1=-1; | |
196 | if(jrLtGeant[ind]!=jrLt[0]) conhwr(jrLtGeant[ind],&incrCount1); | |
197 | jrLtGeant[ind]=-1; | |
198 | } | |
199 | LttcFlagGeant=-1; | |
200 | } | |
201 | ||
202 | ||
203 | //update jrLt and sLt arrays | |
204 | G4int end = 99; | |
205 | if(LttcFlag>=0) | |
206 | end=LttcFlag; | |
207 | ||
208 | // Added by A.Solodkov | |
209 | if (end>=100) { | |
210 | G4cout << "Problems in WrapG1 routine" << G4endl; | |
211 | G4cout << "Index LttcFlag=" << end << " is outside array bounds" << G4endl; | |
212 | G4cout << "Better to stop immediately !" << G4endl; | |
213 | exit(1); | |
214 | } | |
215 | //jrLt re-initialization with -1 (jrLt[0] is already set) | |
216 | for(G4int vv=1;vv<=end;vv++) | |
217 | jrLt[vv]=-1; | |
218 | //sLt re-initialization | |
219 | for(G4int vs=0;vs<=end;vs++) | |
220 | sLt[vs]=0; | |
221 | ||
222 | LttcFlag=0; | |
223 | ||
224 | ||
225 | ///////////////////// COMPUTE STEP //////////////////////// | |
226 | //update Navigator private flags and voxel stack if point is | |
227 | //changed from last located point (otherwise troubles come | |
228 | //when fluka changes location or particle because G4 computes | |
229 | //from last located point). | |
230 | G4ThreeVector newLocalPoint = ptrNavig->ComputeLocalPoint(partLoc); | |
231 | G4double moveLenSq = (newLocalPoint-oldLocalPoint).mag2(); | |
232 | if(moveLenSq>=kCarTolerance*kCarTolerance) | |
233 | ptrNavig->LocateGlobalPointWithinVolume(partLoc); | |
234 | ||
235 | ||
236 | //compute step and new location | |
237 | newReg = oldReg; | |
238 | G4double appStep = 0; | |
239 | G4double safety = 0; | |
240 | G4bool onBoundaryRound = false; | |
241 | G4bool crossBound = false; | |
242 | G4double physStepTmp = physStep; | |
243 | G4bool flagError = false; | |
244 | ||
245 | while( (oldReg==newReg && appStep<physStepTmp) || onBoundaryRound ) { | |
246 | #ifdef G4GEOMETRY_DEBUG | |
247 | G4cout << "* Computing step..." << G4endl; | |
248 | #endif | |
249 | ||
250 | //update variables | |
251 | oldReg = newReg; | |
252 | ||
253 | if(onBoundaryRound) { | |
254 | physStepTmp=10.e-10; | |
255 | //compute step and location: returns newReg | |
256 | newReg=StepAndLocation(partLoc,pVec,physStepTmp,appStep, | |
257 | safety,onBoundaryRound,flagError,oldReg); | |
258 | physStepTmp=0.; | |
259 | crossBound=true; | |
260 | } | |
261 | else { | |
262 | physStepTmp -= appStep; | |
263 | //compute step and location: returns newReg | |
264 | newReg=StepAndLocation(partLoc,pVec,physStepTmp,appStep, | |
265 | safety,onBoundaryRound,flagError,oldReg); | |
266 | } | |
267 | ||
268 | G4bool EqualHist; | |
269 | EqualHist = EqualHistories(ptrTouchHist->GetHistory(), | |
270 | ptrGeoInit->GetTempNavHist()->GetHistory()); | |
271 | ||
272 | if(!EqualHist && flagError) { | |
273 | pVec=-pVec; | |
274 | newReg=StepAndLocation(partLoc,pVec,physStepTmp,appStep, | |
275 | safety,onBoundaryRound,flagError,oldReg); | |
276 | pVec=-pVec; | |
277 | physStepTmp+=1; | |
278 | newReg=StepAndLocation(partLoc,pVec,physStepTmp,appStep, | |
279 | safety,onBoundaryRound,flagError,oldReg); | |
280 | ||
281 | EqualHist = EqualHistories(ptrTouchHist->GetHistory(), | |
282 | ptrGeoInit->GetTempNavHist()->GetHistory()); | |
283 | } | |
284 | ||
285 | //update sLt | |
286 | G4double pas = G4double(appStep); | |
287 | sLt[LttcFlag] += pas/cm; | |
288 | safety = (oldReg!=newReg)?0.:safety; | |
289 | ||
290 | if(!EqualHist) { | |
291 | //end-step point is on boundary between volumes; | |
292 | //==> update step-histories, save new NavHistWithCount | |
293 | //and save its pointer in jrLt; save step in sLt. | |
294 | ||
295 | //set onBoundaryRound=false to avoid re-compute step! | |
296 | onBoundaryRound=false; | |
297 | ||
298 | #ifdef G4GEOMETRY_DEBUG | |
299 | G4cout << "* History is changed!" << G4endl; | |
300 | G4cout << "* updating step-histories, jrLt, LttcFlag" << G4endl; | |
301 | #endif | |
302 | ||
303 | ptrGeoInit->UpdateHistories(ptrTouchHist->GetHistory(),2); | |
304 | ||
305 | LttcFlag += 1; | |
306 | ||
307 | #ifdef G4GEOMETRY_DEBUG | |
308 | G4cout << "* ISVHWR call to store new NavHistWithCount in jrLt" | |
309 | << G4endl; | |
310 | #endif | |
311 | ||
312 | jrLt[LttcFlag] = isvhwr(0,0); | |
313 | ||
314 | #ifdef G4GEOMETRY_DEBUG | |
315 | G4cout << "* CONHWR call to increment counter" << G4endl; | |
316 | #endif | |
317 | G4int incrCount3=1; | |
318 | conhwr(jrLt[LttcFlag],&incrCount3); | |
319 | ||
320 | sLt[LttcFlag] = sLt[LttcFlag-1]; | |
321 | } | |
322 | } | |
323 | ||
324 | ////////////////////// OUTPUT ////////////////////////// | |
325 | //If back-scattering occured, and fluka is in the wrong region, return -3. | |
326 | //(N. B. Boundary between replicans are not seen when physStep=distance | |
327 | //particle-boundary: in this case step=kInfinity and history is unchanged. | |
328 | //Following step is =0 and then history changes.) | |
329 | if(nascFlag<0 && !appStep && physStep && newReg!=oldReg && !crossBound) { | |
330 | //don't need to compare histories because boundary between | |
331 | //identical volumes in different replicans are not seen | |
332 | #ifdef G4GEOMETRY_DEBUG | |
333 | G4cout << "* Back-scattering!" << G4endl; | |
334 | #endif | |
335 | ||
336 | newReg=-3; | |
337 | } | |
338 | else if(newRegErr<0) | |
339 | newReg=newRegErr; | |
340 | ||
341 | ||
342 | //compute output variables (in cm.!) | |
343 | //final step | |
344 | retStep = sLt[LttcFlag]; | |
345 | ||
346 | //safety (Fluka sottracts a bit to safety to be sure | |
347 | //not to jump on a boundary) | |
348 | G4double s = G4double(safety); | |
349 | s -= s*3.0e-09; | |
350 | saf=s/cm; | |
351 | ||
352 | //update wrapper utility variables | |
353 | //copy jrLt in jrLtGeant | |
354 | G4int start=0; | |
355 | if(haveHistNb!=-1 && LttcFlagGeant!=-1) | |
356 | start=1; | |
357 | for(G4int lt=start;lt<=LttcFlag;lt++) | |
358 | jrLtGeant[LttcFlagGeant+1+lt-start]=jrLt[lt]; | |
359 | LttcFlagGeant+=(1+LttcFlag-start); | |
360 | newLttc = jrLt[LttcFlag]; | |
361 | ptrGeoInit->SetLttcFlagGeant(LttcFlagGeant); | |
362 | ||
363 | partLocOld=partLoc; | |
364 | oldLocalPoint = ptrNavig->ComputeLocalPoint(partLocOld); | |
365 | ||
366 | //compute new position | |
367 | G4ThreeVector oldPos = G4ThreeVector(pSx,pSy,pSz); | |
368 | G4ThreeVector newPos = oldPos + retStep*pVec; | |
369 | ||
370 | #ifdef G4GEOMETRY_DEBUG | |
371 | G4cout << "New position: " << newPos << G4endl; | |
372 | G4cout << "Output region: " << newReg << G4endl; | |
373 | G4cout << "G4 safety (cm): " << (safety*0.1) << G4endl; | |
374 | G4cout << "Fluka safety (cm): " << saf << G4endl; | |
375 | G4cout << "Approved step: " << retStep << G4endl; | |
376 | G4cout << "LttcFlag = " << LttcFlag << G4endl; | |
377 | for(G4int i=0;i<=LttcFlag+1;i++) { | |
378 | G4cout << "jrLt[" << i <<"]=" << jrLt[i] << G4endl; | |
379 | G4cout << "sLt[" << i <<"]=" << sLt[i] << G4endl; | |
380 | } | |
381 | ||
382 | G4cout << "LttcFlagGeant =" << LttcFlagGeant << G4endl; | |
383 | for(G4int ib=0;ib<=LttcFlagGeant+1;ib++) { | |
384 | G4cout << "jrLtGeant[" << ib <<"]=" << jrLtGeant[ib] << G4endl; | |
385 | } | |
386 | G4cout << "newLttc=" << newLttc << G4endl; | |
387 | #endif | |
388 | } | |
389 | ||
390 | ||
391 | ||
392 |