Using TGeo to retrieve the mean material budget between two points (M.Ivanov)
[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
68 // change because of geant4 6.0
69 //static G4ThreeVector oldLocalPoint = 
70 //  ptrNavig->ComputeLocalPoint(partLocOld);
71  static G4ThreeVector oldLocalPoint =
72     ptrNavig->GetGlobalToLocalTransform().TransformPoint(partLocOld);
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).
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
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;
351 //======================= Endre print ===================  
352 printf("\nWrapG1 mreg=%d",newReg);
353 //======================= Endre print ===================  
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;
377
378 // change because of geant4 6.0
379 //oldLocalPoint = ptrNavig->ComputeLocalPoint(partLocOld);
380   oldLocalPoint =
381     ptrNavig->GetGlobalToLocalTransform().TransformPoint(partLocOld);
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