First implementation. Needs cleanup.
[u/mrichter/AliRoot.git] / Flugg / WrapUtils.cxx
1 #include "WrapUtils.hh"
2
3 #include "FGeometryInit.hh"
4 #include <iostream.h>
5 #include <iomanip.h>
6
7 ////////////////////////////////////////////////////////////////////////
8 // StepAndLocation
9 ////////////////////////////////////////////////////////////////////////
10
11 G4int StepAndLocation(G4ThreeVector &partLoc, const G4ThreeVector &pVec,
12                       const G4double &physStep, G4double &appStep,
13                       G4double &safety, G4bool & onBound, 
14                       G4bool & flagErr, const G4int & oldReg)
15 {
16   //NB onBound=true in the particular case when particle is in boundary 
17   //rounding BUT geometry hasn't crossed it.
18   //flagErr=true when particle is on boundary but step=infinity, 
19   //newSafety>10.e-10 and lLocate function locates end step point in
20   // a new region. This is absurb!
21   
22   
23   //Geoinit, Navigator, TouchableHistory, etc. pointers
24   static FGeometryInit * ptrGeoInit = FGeometryInit::GetInstance();
25   FluggNavigator * ptrNavig = ptrGeoInit->getNavigatorForTracking();
26   G4TouchableHistory * ptrTouchableHistory = ptrGeoInit->GetTouchableHistory();
27   G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
28   
29   //compute step
30   G4double step = ptrNavig->ComputeStep(partLoc,pVec,physStep,safety);
31   
32   //compute approved step
33   if(step<physStep) {
34     //NB:to check if point is really on boundary: 
35     //see G4AuxiliaryNavServices.hh
36 #ifdef G4GEOMETRY_DEBUG
37     G4cout << "* Boundary crossing!" << G4endl; 
38 #endif
39     
40     appStep=step;
41     ptrNavig->SetGeometricallyLimitedStep();
42     
43     //if SetGeometricallyLimitedStep() is called, next 
44     //LocateGlobalPointAndUpdateTouchable(...,true) will 
45     //start searching from the end of the previous step, 
46     //otherwise will start from last located point. 
47     //(e.g.PropagatorInField and SteppingManger.
48   }
49   else {
50     //step=kInfinity (the nearest boundary is >physStep)
51     //Fluka always wants step lenght approved, so:
52     appStep=physStep;
53   }
54   
55   //new position
56   partLoc+=(pVec*appStep);
57   
58   //locate point and update touchable history
59   ptrNavig->LocateGlobalPointAndUpdateTouchable(partLoc,0,
60                                                 ptrTouchableHistory,true);
61   
62   G4VPhysicalVolume * touchvolume=ptrTouchableHistory->GetVolume();
63   
64   //if volume not found, out of mother volume: 
65   //returns [number of volumes]+1
66   G4int newReg;
67   if(!touchvolume) {
68 #ifdef G4GEOMETRY_DEBUG
69     G4cout << "Out of mother volume! (G1WR)"  << G4endl;
70 #endif
71     G4int numVol = G4int(pVolStore->size());
72     newReg = numVol + 1;
73   }
74   else {
75 #ifdef G4GEOMETRY_DEBUG
76     G4cout <<"Volume after step: " <<touchvolume->GetName() << G4endl;
77 #endif
78     
79     //compute new volume index
80     G4int volIndex = ~0;
81     for (unsigned int i=0; i<pVolStore->size(); i++)
82       if ((*pVolStore)[i] == touchvolume) 
83         volIndex = i;
84     if (volIndex==(~0)) {
85       G4cerr << "ERROR in FLUGG: Problem in routine WrapG1 trying to find volume after step" << G4endl;
86       exit(-999);
87     }
88     newReg=volIndex+1;   
89   }
90   
91   onBound=false;
92   flagErr=false;
93   //check if position is in a boundary rounding while volume isn't changed
94   if(step==kInfinity || step==physStep) {
95     G4ThreeVector globalpoint = ptrNavig->GetLocalToGlobalTransform().
96       TransformPoint(ptrNavig->GetCurrentLocalCoordinate());
97     
98     //compute new safety
99     G4double newSafetydbg = ptrNavig->ComputeSafety(globalpoint,DBL_MAX );
100     
101 #ifdef G4GEOMETRY_DEBUG
102     G4cout << "|From StepAndLocation:"  << G4endl;
103     G4cout << "|step=" << step << G4endl;
104     G4cout << "|phystep=" << physStep << G4endl;
105     G4cout << "|safety=" << safety << G4endl;
106     G4cout << "|newSafetydbg =" << newSafetydbg << G4endl;
107 #endif
108     
109     if(newSafetydbg<1.0e-10) 
110       onBound=true;
111     else if(newReg != oldReg) { 
112 #ifdef G4GEOMETRY_DEBUG
113       G4cout << "New Volume but ComputeStep didn't notice!" << G4endl;
114 #endif
115       flagErr=true;
116     } 
117   }
118   
119   //return
120   return newReg;
121 }
122
123
124
125
126 ////////////////////////////////////////////////////////////////////////
127 // EqualHistories
128 ////////////////////////////////////////////////////////////////////////
129
130 bool EqualHistories(const G4NavigationHistory* ptrFirstHist,
131                     const G4NavigationHistory* ptrSecHist)
132 {
133 #ifdef G4GEOMETRY_DEBUG
134   G4cout << "* Testing if histories are equal..."  << G4endl;
135 #endif
136      
137   G4int depth1 = ptrFirstHist->GetDepth();
138   G4int depth2 = ptrSecHist->GetDepth();
139   
140   if(depth1!=depth2) 
141     return false;
142   
143   for(G4int w=0;w<=depth1;w++) {
144     if (*(ptrFirstHist->GetVolume(w))==*(ptrSecHist->GetVolume(w))) {
145       //kNormal volume
146       if(ptrFirstHist->GetVolumeType(w) == kNormal) {
147         if ( ptrFirstHist->GetVolume(w)->GetCopyNo() !=
148              ptrSecHist->GetVolume(w)->GetCopyNo() )
149           return false;
150       }
151       
152       //Replica or Parametric volume
153       else {
154         if ( ptrFirstHist->GetReplicaNo(w) !=
155              ptrSecHist->GetReplicaNo(w) )
156           return false;
157       }
158     }
159     else
160       return false;     
161   }
162   
163 #ifdef G4GEOMETRY_DEBUG
164   G4cout << "Histories are equal!" << G4endl;   
165 #endif
166   
167   return true;
168 }
169
170
171
172 ////////////////////////////////////////////////////////////////////////
173 // PrintHeader
174 ////////////////////////////////////////////////////////////////////////
175 ostream& PrintHeader(ostream& os, const char* title) {
176   os << "*\n" << "*\n" << "*\n";
177   os << "*********************  " << title << " *********************\n"
178      << "*\n";
179   os << "*...+....1....+....2....+....3....+....4....+....5....+....6....+....7..."
180      << endl;
181   os << "*" << endl;
182
183   return os;
184 }
185
186 ////////////////////////////////////////////////////////////////////////
187 // PrintMaterial
188 ////////////////////////////////////////////////////////////////////////
189 ostream& PrintMaterial(ostream& os, const char* title,
190                        G4double Z, G4double A,
191                        G4double density,
192                        G4double index,
193                        G4double N,
194                        const char* name) {
195
196   os << setw10 << title;
197
198   os.setf(0,G4std::ios::floatfield);
199   if (Z < 0)
200     os << setw10 << " ";
201   else
202     os << setw10 
203        << setfixed
204        << G4std::setprecision(1) 
205        << Z;
206   
207   if (A < 0)
208     os << setw10 << " ";
209   else
210     os << setw10 << G4std::setprecision(3)
211        << A;
212
213   os.setf(0,G4std::ios::floatfield);
214   os << setw10 
215      << setscientific
216      << G4std::setprecision(3) 
217      << density;
218
219   os.setf(0,G4std::ios::floatfield);
220   os << setw10 
221      << setfixed
222      << G4std::setprecision(1) 
223      << index;
224
225   os << setw10 << " ";
226
227   if (N < 0)
228     os << setw10 << " ";
229   else
230     os << setw10 << N;
231
232   os << name << endl;
233
234   return os;
235 }
236
237
238 ////////////////////////////////////////////////////////////////////////
239 // PrintCompund
240 ////////////////////////////////////////////////////////////////////////
241 ostream& PrintCompound(ostream& os, const char* title,
242                        G4int count,
243                        const char* name,
244                        G4double fraction,
245                        G4double index) {
246
247   
248   if(count==3) {
249     os << name << endl;
250     os << setw10 << "COMPOUND  ";
251   }
252   
253   os.setf(0,G4std::ios::floatfield);
254   os << setw10
255      << setscientific
256      << G4std::setprecision(2)
257      << fraction;
258   os.setf(0,G4std::ios::floatfield);
259   os << setw10
260      << setfixed
261      << G4std::setprecision(1)
262      << index;
263
264   return os;
265 }
266