Including St2 in the new geometry segmentation (Christian)
[u/mrichter/AliRoot.git] / Flugg / WrapG1.cxx
CommitLineData
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
41void 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 ===================
352printf("\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