]>
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; | |
0edf14e7 | 67 | |
68 | // change because of geant4 6.0 | |
69 | //static G4ThreeVector oldLocalPoint = | |
70 | // ptrNavig->ComputeLocalPoint(partLocOld); | |
71 | static G4ThreeVector oldLocalPoint = | |
72 | ptrNavig->GetGlobalToLocalTransform().TransformPoint(partLocOld); | |
26911512 | 73 | |
74 | G4ThreeVector pVec(pV[0],pV[1],pV[2]); | |
75 | const G4double physStep=G4double(propStep*10.); | |
76 | ||
77 | G4int LttcFlagGeant = ptrGeoInit->GetLttcFlagGeant(); | |
78 | ||
79 | ||
80 | #ifdef G4GEOMETRY_DEBUG | |
81 | G4cout.precision(10); | |
82 | G4cout << "Position (cm):" << pSx << "," << pSy << "," << pSz << G4endl; | |
83 | G4cout << "Direction: " << pVec << G4endl; | |
84 | G4cout << "Proposed step :"<< propStep << G4endl; | |
85 | #endif | |
86 | ||
87 | ||
88 | ||
89 | ///////////////////// FLUKA/G4 REGIONS COMPARISON ////////////////////// | |
90 | //get oldReg pointer and G4volume of previous step end-point pointer | |
91 | G4VPhysicalVolume * ptrOldReg = (*pVolStore)[oldReg-1]; | |
92 | G4VPhysicalVolume * lastVolume = ptrTouchHist->GetVolume(); | |
93 | G4int touVolCpNb = 0; | |
94 | if (lastVolume) | |
95 | touVolCpNb = lastVolume->GetCopyNo(); | |
96 | ||
97 | #ifdef G4GEOMETRY_DEBUG | |
98 | G4int fluVolCpNb = ptrOldReg->GetCopyNo(); | |
99 | G4cout << "Fluka volume before step: " << ptrOldReg->GetName() | |
100 | << "," << fluVolCpNb<< G4endl; | |
101 | if(lastVolume) | |
102 | G4cout << "G4 Touch Hist volume: " | |
103 | << lastVolume->GetName() << "," << touVolCpNb << G4endl; | |
104 | G4cout << "------------------------------------------------" << G4endl; | |
105 | #endif | |
106 | ||
107 | ||
108 | //if volume is changed, this is a new particle tracking, or fluka tries | |
109 | //to reach a boundary more softly with the same particle. In this case | |
110 | //fluka restart tracking from old history, in general. For tracking in | |
111 | //lattice-volume fluka goes back to one of the previous lattice volumes. | |
112 | //Check if ptrOldReg is equal to old history, or to one of the N lattice | |
113 | //volumes stored in jrLt. Then reinitialise step-histories! Otherwise | |
114 | //must relocate. | |
115 | //NB. jrLtGeant stores lattice volume histories until LttcFlag==-1, | |
116 | //then all histories are checked and deleted, and array is reinitialised | |
117 | ||
118 | ||
119 | G4int haveHistNb = -1; | |
120 | G4int newRegErr=0; | |
121 | G4int indHist = LttcFlagGeant; | |
122 | G4int * jrLtGeant = ptrGeoInit->GetJrLtGeantArray(); | |
123 | ||
124 | ||
125 | G4NavigationHistory* ptrLttcHist=0; | |
126 | if(oldLttc) | |
127 | ptrLttcHist = reinterpret_cast<NavHistWithCount*>(oldLttc)->GetNavHistPtr(); | |
128 | ||
129 | ||
130 | while(indHist>=0 && haveHistNb==-1) { | |
131 | #ifdef G4GEOMETRY_DEBUG | |
132 | G4cout<<"Searching in jrLtc...."<<G4endl; | |
133 | #endif | |
134 | if(oldLttc==jrLtGeant[indHist]) | |
135 | haveHistNb=indHist; | |
136 | indHist-=1; | |
137 | } | |
138 | ||
139 | if(haveHistNb!=-1) { | |
140 | //fluka history found in jrLtGeant | |
141 | if(haveHistNb<LttcFlagGeant) { | |
142 | #ifdef G4GEOMETRY_DEBUG | |
143 | G4cout<<"* Fluka reaches boundary more softly..."<<G4endl; | |
144 | G4cout<<"* Re-initializing FluggNavigator history"<<G4endl; | |
145 | G4cout<<"* and updating step-histories"<<G4endl; | |
146 | #endif | |
147 | ||
148 | ptrNavig->UpdateNavigatorHistory(ptrLttcHist); | |
149 | ptrGeoInit->UpdateHistories(ptrLttcHist,0); | |
150 | } | |
151 | #ifdef G4GEOMETRY_DEBUG | |
152 | if(haveHistNb==LttcFlagGeant) | |
153 | G4cout << "Continuing step...." << G4endl; | |
154 | #endif | |
155 | jrLt[0]=oldLttc; | |
156 | } | |
157 | else { | |
158 | //not found fluka history in jrLttGeant! | |
159 | G4cout << "* ERROR! Geometry not correctly initialised in fluka history!" | |
160 | << G4endl; | |
161 | ||
162 | //relocation! | |
163 | ptrNavig->LocateGlobalPointAndUpdateTouchable(partLoc,0,ptrTouchHist,true); | |
164 | ||
165 | G4cout << "* ATTENTION: point relocation in: " | |
166 | << ptrTouchHist->GetVolume()->GetName() << G4endl; | |
167 | ||
168 | ptrGeoInit->UpdateHistories(ptrTouchHist->GetHistory(),1); | |
169 | ||
170 | if(ptrTouchHist->GetVolume() != ptrOldReg) { | |
171 | G4cout << "* ERROR! Point not in fluka volume!" << G4endl; | |
172 | newRegErr=-3; | |
173 | } | |
174 | ||
175 | //save new history in jrLt[0] and increment its counter | |
176 | #ifdef G4GEOMETRY_DEBUG | |
177 | G4cout << "* ISVHWR call to store new NavHistWithCount in jrLt[0]" | |
178 | << G4endl; | |
179 | #endif | |
180 | jrLt[0]=isvhwr(0,0); | |
181 | ||
182 | #ifdef G4GEOMETRY_DEBUG | |
183 | G4cout << "* CONHWR call to increment counter" << G4endl; | |
184 | #endif | |
185 | G4int incrCount2=1; | |
186 | conhwr(jrLt[0],&incrCount2); | |
187 | } | |
188 | ||
189 | ||
190 | //jrLtGeant - history check: decrement counter and delete | |
191 | //histories, if LttcFlag=-1, then reinitialise array with -1. | |
192 | if(LttcFlag==-1 && LttcFlagGeant>=0) { | |
193 | #ifdef G4GEOMETRY_DEBUG | |
194 | G4cout << "* CONHWR call to check and delete histories in jrLtGeant[]:" | |
195 | << G4endl; | |
196 | #endif | |
197 | ||
198 | for(G4int ind=0; ind<=LttcFlagGeant; ind++) { | |
199 | G4int incrCount1=-1; | |
200 | if(jrLtGeant[ind]!=jrLt[0]) conhwr(jrLtGeant[ind],&incrCount1); | |
201 | jrLtGeant[ind]=-1; | |
202 | } | |
203 | LttcFlagGeant=-1; | |
204 | } | |
205 | ||
206 | ||
207 | //update jrLt and sLt arrays | |
208 | G4int end = 99; | |
209 | if(LttcFlag>=0) | |
210 | end=LttcFlag; | |
211 | ||
212 | // Added by A.Solodkov | |
213 | if (end>=100) { | |
214 | G4cout << "Problems in WrapG1 routine" << G4endl; | |
215 | G4cout << "Index LttcFlag=" << end << " is outside array bounds" << G4endl; | |
216 | G4cout << "Better to stop immediately !" << G4endl; | |
217 | exit(1); | |
218 | } | |
219 | //jrLt re-initialization with -1 (jrLt[0] is already set) | |
220 | for(G4int vv=1;vv<=end;vv++) | |
221 | jrLt[vv]=-1; | |
222 | //sLt re-initialization | |
223 | for(G4int vs=0;vs<=end;vs++) | |
224 | sLt[vs]=0; | |
225 | ||
226 | LttcFlag=0; | |
227 | ||
228 | ||
229 | ///////////////////// COMPUTE STEP //////////////////////// | |
230 | //update Navigator private flags and voxel stack if point is | |
231 | //changed from last located point (otherwise troubles come | |
232 | //when fluka changes location or particle because G4 computes | |
233 | //from last located point). | |
0edf14e7 | 234 | |
235 | // change because of geant4 6.0 | |
236 | //G4ThreeVector newLocalPoint = ptrNavig->ComputeLocalPoint(partLoc); | |
237 | G4ThreeVector newLocalPoint = | |
238 | ptrNavig->GetGlobalToLocalTransform().TransformPoint(partLoc); | |
239 | ||
240 | ||
241 | ||
26911512 | 242 | G4double moveLenSq = (newLocalPoint-oldLocalPoint).mag2(); |
243 | if(moveLenSq>=kCarTolerance*kCarTolerance) | |
244 | ptrNavig->LocateGlobalPointWithinVolume(partLoc); | |
245 | ||
246 | ||
247 | //compute step and new location | |
248 | newReg = oldReg; | |
249 | G4double appStep = 0; | |
250 | G4double safety = 0; | |
251 | G4bool onBoundaryRound = false; | |
252 | G4bool crossBound = false; | |
253 | G4double physStepTmp = physStep; | |
254 | G4bool flagError = false; | |
255 | ||
256 | while( (oldReg==newReg && appStep<physStepTmp) || onBoundaryRound ) { | |
257 | #ifdef G4GEOMETRY_DEBUG | |
258 | G4cout << "* Computing step..." << G4endl; | |
259 | #endif | |
260 | ||
261 | //update variables | |
262 | oldReg = newReg; | |
263 | ||
264 | if(onBoundaryRound) { | |
265 | physStepTmp=10.e-10; | |
266 | //compute step and location: returns newReg | |
267 | newReg=StepAndLocation(partLoc,pVec,physStepTmp,appStep, | |
268 | safety,onBoundaryRound,flagError,oldReg); | |
269 | physStepTmp=0.; | |
270 | crossBound=true; | |
271 | } | |
272 | else { | |
273 | physStepTmp -= appStep; | |
274 | //compute step and location: returns newReg | |
275 | newReg=StepAndLocation(partLoc,pVec,physStepTmp,appStep, | |
276 | safety,onBoundaryRound,flagError,oldReg); | |
277 | } | |
278 | ||
279 | G4bool EqualHist; | |
280 | EqualHist = EqualHistories(ptrTouchHist->GetHistory(), | |
281 | ptrGeoInit->GetTempNavHist()->GetHistory()); | |
282 | ||
283 | if(!EqualHist && flagError) { | |
284 | pVec=-pVec; | |
285 | newReg=StepAndLocation(partLoc,pVec,physStepTmp,appStep, | |
286 | safety,onBoundaryRound,flagError,oldReg); | |
287 | pVec=-pVec; | |
288 | physStepTmp+=1; | |
289 | newReg=StepAndLocation(partLoc,pVec,physStepTmp,appStep, | |
290 | safety,onBoundaryRound,flagError,oldReg); | |
291 | ||
292 | EqualHist = EqualHistories(ptrTouchHist->GetHistory(), | |
293 | ptrGeoInit->GetTempNavHist()->GetHistory()); | |
294 | } | |
295 | ||
296 | //update sLt | |
297 | G4double pas = G4double(appStep); | |
298 | sLt[LttcFlag] += pas/cm; | |
299 | safety = (oldReg!=newReg)?0.:safety; | |
300 | ||
301 | if(!EqualHist) { | |
302 | //end-step point is on boundary between volumes; | |
303 | //==> update step-histories, save new NavHistWithCount | |
304 | //and save its pointer in jrLt; save step in sLt. | |
305 | ||
306 | //set onBoundaryRound=false to avoid re-compute step! | |
307 | onBoundaryRound=false; | |
308 | ||
309 | #ifdef G4GEOMETRY_DEBUG | |
310 | G4cout << "* History is changed!" << G4endl; | |
311 | G4cout << "* updating step-histories, jrLt, LttcFlag" << G4endl; | |
312 | #endif | |
313 | ||
314 | ptrGeoInit->UpdateHistories(ptrTouchHist->GetHistory(),2); | |
315 | ||
316 | LttcFlag += 1; | |
317 | ||
318 | #ifdef G4GEOMETRY_DEBUG | |
319 | G4cout << "* ISVHWR call to store new NavHistWithCount in jrLt" | |
320 | << G4endl; | |
321 | #endif | |
322 | ||
323 | jrLt[LttcFlag] = isvhwr(0,0); | |
324 | ||
325 | #ifdef G4GEOMETRY_DEBUG | |
326 | G4cout << "* CONHWR call to increment counter" << G4endl; | |
327 | #endif | |
328 | G4int incrCount3=1; | |
329 | conhwr(jrLt[LttcFlag],&incrCount3); | |
330 | ||
331 | sLt[LttcFlag] = sLt[LttcFlag-1]; | |
332 | } | |
333 | } | |
334 | ||
335 | ////////////////////// OUTPUT ////////////////////////// | |
336 | //If back-scattering occured, and fluka is in the wrong region, return -3. | |
337 | //(N. B. Boundary between replicans are not seen when physStep=distance | |
338 | //particle-boundary: in this case step=kInfinity and history is unchanged. | |
339 | //Following step is =0 and then history changes.) | |
340 | if(nascFlag<0 && !appStep && physStep && newReg!=oldReg && !crossBound) { | |
341 | //don't need to compare histories because boundary between | |
342 | //identical volumes in different replicans are not seen | |
343 | #ifdef G4GEOMETRY_DEBUG | |
344 | G4cout << "* Back-scattering!" << G4endl; | |
345 | #endif | |
346 | ||
347 | newReg=-3; | |
348 | } | |
349 | else if(newRegErr<0) | |
350 | newReg=newRegErr; | |
0edf14e7 | 351 | //======================= Endre print =================== |
352 | printf("\nWrapG1 mreg=%d",newReg); | |
353 | //======================= Endre print =================== | |
26911512 | 354 | |
355 | //compute output variables (in cm.!) | |
356 | //final step | |
357 | retStep = sLt[LttcFlag]; | |
358 | ||
359 | //safety (Fluka sottracts a bit to safety to be sure | |
360 | //not to jump on a boundary) | |
361 | G4double s = G4double(safety); | |
362 | s -= s*3.0e-09; | |
363 | saf=s/cm; | |
364 | ||
365 | //update wrapper utility variables | |
366 | //copy jrLt in jrLtGeant | |
367 | G4int start=0; | |
368 | if(haveHistNb!=-1 && LttcFlagGeant!=-1) | |
369 | start=1; | |
370 | for(G4int lt=start;lt<=LttcFlag;lt++) | |
371 | jrLtGeant[LttcFlagGeant+1+lt-start]=jrLt[lt]; | |
372 | LttcFlagGeant+=(1+LttcFlag-start); | |
373 | newLttc = jrLt[LttcFlag]; | |
374 | ptrGeoInit->SetLttcFlagGeant(LttcFlagGeant); | |
375 | ||
376 | partLocOld=partLoc; | |
0edf14e7 | 377 | |
378 | // change because of geant4 6.0 | |
379 | //oldLocalPoint = ptrNavig->ComputeLocalPoint(partLocOld); | |
380 | oldLocalPoint = | |
381 | ptrNavig->GetGlobalToLocalTransform().TransformPoint(partLocOld); | |
26911512 | 382 | |
383 | //compute new position | |
384 | G4ThreeVector oldPos = G4ThreeVector(pSx,pSy,pSz); | |
385 | G4ThreeVector newPos = oldPos + retStep*pVec; | |
386 | ||
387 | #ifdef G4GEOMETRY_DEBUG | |
388 | G4cout << "New position: " << newPos << G4endl; | |
389 | G4cout << "Output region: " << newReg << G4endl; | |
390 | G4cout << "G4 safety (cm): " << (safety*0.1) << G4endl; | |
391 | G4cout << "Fluka safety (cm): " << saf << G4endl; | |
392 | G4cout << "Approved step: " << retStep << G4endl; | |
393 | G4cout << "LttcFlag = " << LttcFlag << G4endl; | |
394 | for(G4int i=0;i<=LttcFlag+1;i++) { | |
395 | G4cout << "jrLt[" << i <<"]=" << jrLt[i] << G4endl; | |
396 | G4cout << "sLt[" << i <<"]=" << sLt[i] << G4endl; | |
397 | } | |
398 | ||
399 | G4cout << "LttcFlagGeant =" << LttcFlagGeant << G4endl; | |
400 | for(G4int ib=0;ib<=LttcFlagGeant+1;ib++) { | |
401 | G4cout << "jrLtGeant[" << ib <<"]=" << jrLtGeant[ib] << G4endl; | |
402 | } | |
403 | G4cout << "newLttc=" << newLttc << G4endl; | |
404 | #endif | |
405 | } | |
406 | ||
407 | ||
408 | ||
409 |