4 ///////////////////////////////////////////////////////////////////
6 // WrapG1.hh - Sara Vanini
8 // Wrapper for geometry tracking: returns approved step of
9 // particle and alla variables that fluka G1 computes.
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.
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);
28 /////////////////////////////////////////////////////////////////////
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"
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)
48 #ifdef G4GEOMETRY_DEBUG
49 G4cout<<"============= G1WR =============="<<G4endl;
52 //static G4int count=0;
54 //G4cout<<"contatore G1="<<count<<G4endl;
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();
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);
70 G4ThreeVector pVec(pV[0],pV[1],pV[2]);
71 const G4double physStep=G4double(propStep*10.);
73 G4int LttcFlagGeant = ptrGeoInit->GetLttcFlagGeant();
76 #ifdef G4GEOMETRY_DEBUG
78 G4cout << "Position (cm):" << pSx << "," << pSy << "," << pSz << G4endl;
79 G4cout << "Direction: " << pVec << G4endl;
80 G4cout << "Proposed step :"<< propStep << G4endl;
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();
91 touVolCpNb = lastVolume->GetCopyNo();
93 #ifdef G4GEOMETRY_DEBUG
94 G4int fluVolCpNb = ptrOldReg->GetCopyNo();
95 G4cout << "Fluka volume before step: " << ptrOldReg->GetName()
96 << "," << fluVolCpNb<< G4endl;
98 G4cout << "G4 Touch Hist volume: "
99 << lastVolume->GetName() << "," << touVolCpNb << G4endl;
100 G4cout << "------------------------------------------------" << G4endl;
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
111 //NB. jrLtGeant stores lattice volume histories until LttcFlag==-1,
112 //then all histories are checked and deleted, and array is reinitialised
115 G4int haveHistNb = -1;
117 G4int indHist = LttcFlagGeant;
118 G4int * jrLtGeant = ptrGeoInit->GetJrLtGeantArray();
121 G4NavigationHistory* ptrLttcHist=0;
123 ptrLttcHist = reinterpret_cast<NavHistWithCount*>(oldLttc)->GetNavHistPtr();
126 while(indHist>=0 && haveHistNb==-1) {
127 #ifdef G4GEOMETRY_DEBUG
128 G4cout<<"Searching in jrLtc...."<<G4endl;
130 if(oldLttc==jrLtGeant[indHist])
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;
144 ptrNavig->UpdateNavigatorHistory(ptrLttcHist);
145 ptrGeoInit->UpdateHistories(ptrLttcHist,0);
147 #ifdef G4GEOMETRY_DEBUG
148 if(haveHistNb==LttcFlagGeant)
149 G4cout << "Continuing step...." << G4endl;
154 //not found fluka history in jrLttGeant!
155 G4cout << "* ERROR! Geometry not correctly initialised in fluka history!"
159 ptrNavig->LocateGlobalPointAndUpdateTouchable(partLoc,0,ptrTouchHist,true);
161 G4cout << "* ATTENTION: point relocation in: "
162 << ptrTouchHist->GetVolume()->GetName() << G4endl;
164 ptrGeoInit->UpdateHistories(ptrTouchHist->GetHistory(),1);
166 if(ptrTouchHist->GetVolume() != ptrOldReg) {
167 G4cout << "* ERROR! Point not in fluka volume!" << G4endl;
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]"
178 #ifdef G4GEOMETRY_DEBUG
179 G4cout << "* CONHWR call to increment counter" << G4endl;
182 conhwr(jrLt[0],&incrCount2);
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[]:"
194 for(G4int ind=0; ind<=LttcFlagGeant; ind++) {
196 if(jrLtGeant[ind]!=jrLt[0]) conhwr(jrLtGeant[ind],&incrCount1);
203 //update jrLt and sLt arrays
208 // Added by A.Solodkov
210 G4cout << "Problems in WrapG1 routine" << G4endl;
211 G4cout << "Index LttcFlag=" << end << " is outside array bounds" << G4endl;
212 G4cout << "Better to stop immediately !" << G4endl;
215 //jrLt re-initialization with -1 (jrLt[0] is already set)
216 for(G4int vv=1;vv<=end;vv++)
218 //sLt re-initialization
219 for(G4int vs=0;vs<=end;vs++)
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);
236 //compute step and new location
238 G4double appStep = 0;
240 G4bool onBoundaryRound = false;
241 G4bool crossBound = false;
242 G4double physStepTmp = physStep;
243 G4bool flagError = false;
245 while( (oldReg==newReg && appStep<physStepTmp) || onBoundaryRound ) {
246 #ifdef G4GEOMETRY_DEBUG
247 G4cout << "* Computing step..." << G4endl;
253 if(onBoundaryRound) {
255 //compute step and location: returns newReg
256 newReg=StepAndLocation(partLoc,pVec,physStepTmp,appStep,
257 safety,onBoundaryRound,flagError,oldReg);
262 physStepTmp -= appStep;
263 //compute step and location: returns newReg
264 newReg=StepAndLocation(partLoc,pVec,physStepTmp,appStep,
265 safety,onBoundaryRound,flagError,oldReg);
269 EqualHist = EqualHistories(ptrTouchHist->GetHistory(),
270 ptrGeoInit->GetTempNavHist()->GetHistory());
272 if(!EqualHist && flagError) {
274 newReg=StepAndLocation(partLoc,pVec,physStepTmp,appStep,
275 safety,onBoundaryRound,flagError,oldReg);
278 newReg=StepAndLocation(partLoc,pVec,physStepTmp,appStep,
279 safety,onBoundaryRound,flagError,oldReg);
281 EqualHist = EqualHistories(ptrTouchHist->GetHistory(),
282 ptrGeoInit->GetTempNavHist()->GetHistory());
286 G4double pas = G4double(appStep);
287 sLt[LttcFlag] += pas/cm;
288 safety = (oldReg!=newReg)?0.:safety;
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.
295 //set onBoundaryRound=false to avoid re-compute step!
296 onBoundaryRound=false;
298 #ifdef G4GEOMETRY_DEBUG
299 G4cout << "* History is changed!" << G4endl;
300 G4cout << "* updating step-histories, jrLt, LttcFlag" << G4endl;
303 ptrGeoInit->UpdateHistories(ptrTouchHist->GetHistory(),2);
307 #ifdef G4GEOMETRY_DEBUG
308 G4cout << "* ISVHWR call to store new NavHistWithCount in jrLt"
312 jrLt[LttcFlag] = isvhwr(0,0);
314 #ifdef G4GEOMETRY_DEBUG
315 G4cout << "* CONHWR call to increment counter" << G4endl;
318 conhwr(jrLt[LttcFlag],&incrCount3);
320 sLt[LttcFlag] = sLt[LttcFlag-1];
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;
342 //compute output variables (in cm.!)
344 retStep = sLt[LttcFlag];
346 //safety (Fluka sottracts a bit to safety to be sure
347 //not to jump on a boundary)
348 G4double s = G4double(safety);
352 //update wrapper utility variables
353 //copy jrLt in jrLtGeant
355 if(haveHistNb!=-1 && LttcFlagGeant!=-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);
364 oldLocalPoint = ptrNavig->ComputeLocalPoint(partLocOld);
366 //compute new position
367 G4ThreeVector oldPos = G4ThreeVector(pSx,pSy,pSz);
368 G4ThreeVector newPos = oldPos + retStep*pVec;
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;
382 G4cout << "LttcFlagGeant =" << LttcFlagGeant << G4endl;
383 for(G4int ib=0;ib<=LttcFlagGeant+1;ib++) {
384 G4cout << "jrLtGeant[" << ib <<"]=" << jrLtGeant[ib] << G4endl;
386 G4cout << "newLttc=" << newLttc << G4endl;