First implementation. Needs cleanup.
[u/mrichter/AliRoot.git] / Flugg / WrapG1.cxx
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