Including St2 in the new geometry segmentation (Christian)
[u/mrichter/AliRoot.git] / Flugg / WrapUtils.cxx
CommitLineData
26911512 1#include "WrapUtils.hh"
26911512 2#include "FGeometryInit.hh"
26911512 3
4////////////////////////////////////////////////////////////////////////
5// StepAndLocation
6////////////////////////////////////////////////////////////////////////
7
8G4int StepAndLocation(G4ThreeVector &partLoc, const G4ThreeVector &pVec,
9 const G4double &physStep, G4double &appStep,
10 G4double &safety, G4bool & onBound,
11 G4bool & flagErr, const G4int & oldReg)
12{
13 //NB onBound=true in the particular case when particle is in boundary
14 //rounding BUT geometry hasn't crossed it.
15 //flagErr=true when particle is on boundary but step=infinity,
16 //newSafety>10.e-10 and lLocate function locates end step point in
17 // a new region. This is absurb!
18
19
20 //Geoinit, Navigator, TouchableHistory, etc. pointers
21 static FGeometryInit * ptrGeoInit = FGeometryInit::GetInstance();
22 FluggNavigator * ptrNavig = ptrGeoInit->getNavigatorForTracking();
23 G4TouchableHistory * ptrTouchableHistory = ptrGeoInit->GetTouchableHistory();
24 G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
25
26 //compute step
27 G4double step = ptrNavig->ComputeStep(partLoc,pVec,physStep,safety);
28
29 //compute approved step
30 if(step<physStep) {
31 //NB:to check if point is really on boundary:
32 //see G4AuxiliaryNavServices.hh
33#ifdef G4GEOMETRY_DEBUG
34 G4cout << "* Boundary crossing!" << G4endl;
35#endif
36
37 appStep=step;
38 ptrNavig->SetGeometricallyLimitedStep();
39
40 //if SetGeometricallyLimitedStep() is called, next
41 //LocateGlobalPointAndUpdateTouchable(...,true) will
42 //start searching from the end of the previous step,
43 //otherwise will start from last located point.
44 //(e.g.PropagatorInField and SteppingManger.
45 }
46 else {
47 //step=kInfinity (the nearest boundary is >physStep)
48 //Fluka always wants step lenght approved, so:
49 appStep=physStep;
50 }
51
52 //new position
53 partLoc+=(pVec*appStep);
54
55 //locate point and update touchable history
56 ptrNavig->LocateGlobalPointAndUpdateTouchable(partLoc,0,
57 ptrTouchableHistory,true);
58
59 G4VPhysicalVolume * touchvolume=ptrTouchableHistory->GetVolume();
60
61 //if volume not found, out of mother volume:
62 //returns [number of volumes]+1
63 G4int newReg;
64 if(!touchvolume) {
65#ifdef G4GEOMETRY_DEBUG
66 G4cout << "Out of mother volume! (G1WR)" << G4endl;
67#endif
68 G4int numVol = G4int(pVolStore->size());
69 newReg = numVol + 1;
70 }
71 else {
72#ifdef G4GEOMETRY_DEBUG
73 G4cout <<"Volume after step: " <<touchvolume->GetName() << G4endl;
74#endif
75
76 //compute new volume index
77 G4int volIndex = ~0;
78 for (unsigned int i=0; i<pVolStore->size(); i++)
79 if ((*pVolStore)[i] == touchvolume)
80 volIndex = i;
81 if (volIndex==(~0)) {
82 G4cerr << "ERROR in FLUGG: Problem in routine WrapG1 trying to find volume after step" << G4endl;
83 exit(-999);
84 }
85 newReg=volIndex+1;
86 }
87
88 onBound=false;
89 flagErr=false;
90 //check if position is in a boundary rounding while volume isn't changed
91 if(step==kInfinity || step==physStep) {
92 G4ThreeVector globalpoint = ptrNavig->GetLocalToGlobalTransform().
93 TransformPoint(ptrNavig->GetCurrentLocalCoordinate());
94
95 //compute new safety
96 G4double newSafetydbg = ptrNavig->ComputeSafety(globalpoint,DBL_MAX );
97
98#ifdef G4GEOMETRY_DEBUG
99 G4cout << "|From StepAndLocation:" << G4endl;
100 G4cout << "|step=" << step << G4endl;
101 G4cout << "|phystep=" << physStep << G4endl;
102 G4cout << "|safety=" << safety << G4endl;
103 G4cout << "|newSafetydbg =" << newSafetydbg << G4endl;
104#endif
105
106 if(newSafetydbg<1.0e-10)
107 onBound=true;
108 else if(newReg != oldReg) {
109#ifdef G4GEOMETRY_DEBUG
110 G4cout << "New Volume but ComputeStep didn't notice!" << G4endl;
111#endif
112 flagErr=true;
113 }
114 }
115
116 //return
117 return newReg;
118}
119
120
121
122
123////////////////////////////////////////////////////////////////////////
124// EqualHistories
125////////////////////////////////////////////////////////////////////////
126
127bool EqualHistories(const G4NavigationHistory* ptrFirstHist,
128 const G4NavigationHistory* ptrSecHist)
129{
130#ifdef G4GEOMETRY_DEBUG
131 G4cout << "* Testing if histories are equal..." << G4endl;
132#endif
133
134 G4int depth1 = ptrFirstHist->GetDepth();
135 G4int depth2 = ptrSecHist->GetDepth();
136
137 if(depth1!=depth2)
138 return false;
139
140 for(G4int w=0;w<=depth1;w++) {
141 if (*(ptrFirstHist->GetVolume(w))==*(ptrSecHist->GetVolume(w))) {
142 //kNormal volume
143 if(ptrFirstHist->GetVolumeType(w) == kNormal) {
144 if ( ptrFirstHist->GetVolume(w)->GetCopyNo() !=
145 ptrSecHist->GetVolume(w)->GetCopyNo() )
146 return false;
147 }
148
149 //Replica or Parametric volume
150 else {
151 if ( ptrFirstHist->GetReplicaNo(w) !=
152 ptrSecHist->GetReplicaNo(w) )
153 return false;
154 }
155 }
156 else
157 return false;
158 }
159
160#ifdef G4GEOMETRY_DEBUG
161 G4cout << "Histories are equal!" << G4endl;
162#endif
163
164 return true;
165}
166
167
168
169////////////////////////////////////////////////////////////////////////
170// PrintHeader
171////////////////////////////////////////////////////////////////////////
0edf14e7 172std::ostream& PrintHeader(std::ostream& os, const char* title) {
26911512 173 os << "*\n" << "*\n" << "*\n";
174 os << "********************* " << title << " *********************\n"
175 << "*\n";
176 os << "*...+....1....+....2....+....3....+....4....+....5....+....6....+....7..."
36c70081 177 << G4endl;
178 os << "*" << G4endl;
26911512 179
180 return os;
181}
182
183////////////////////////////////////////////////////////////////////////
26d97e06 184// ToFlukaString converts an string to something usefull in FLUKA:
185// * Capital letters
186// * Only 8 letters
187// * Replace ' ' by '_'
26911512 188////////////////////////////////////////////////////////////////////////
26d97e06 189G4String ToFlukaString(G4String str) {
190 str = str.substr(0,8);
191 str.toUpper();
192 for (unsigned int i = 0; i < str.length(); i++) {
193 if (str.at(i) == ' ')
194 str.replace(i,1,"_",1);
26911512 195 }
26d97e06 196 return str;
26911512 197}