24 #include <TGeoMedium.h>
25 #include <TGeoMaterial.h>
26 #include <TGeoManager.h>
39 double dirX,
double dirY,
double dirZ){
41 debugOut <<
"TGeoMaterialInterface::initTrack. \n";
42 debugOut <<
"Pos "; TVector3(posX, posY, posZ).Print();
43 debugOut <<
"Dir "; TVector3(dirX, dirY, dirZ).Print();
47 bool result = !gGeoManager->IsSameLocation(posX, posY, posZ, kTRUE);
49 gGeoManager->SetCurrentDirection(dirX, dirY, dirZ);
52 debugOut <<
" TGeoMaterialInterface::initTrack at \n";
53 debugOut <<
" position: "; TVector3(posX, posY, posZ).Print();
54 debugOut <<
" direction: "; TVector3(dirX, dirY, dirZ).Print();
63 TGeoMaterial* mat = gGeoManager->GetCurrentVolume()->GetMedium()->GetMaterial();
71 const M1x7& stateOrig,
76 const double epsilon(1.
E-1);
80 M1x7 state7, oldState7;
81 oldState7 = stateOrig;
83 int stepSign(sMax < 0 ? -1 : 1);
87 const unsigned maxIt = 300;
91 gGeoManager->FindNextBoundary(fabs(sMax) - s);
92 double safety = gGeoManager->GetSafeDistance();
93 double slDist = gGeoManager->GetStep();
101 double step = slDist;
104 debugOut <<
" safety = " << safety <<
"; step, slDist = " << step <<
"\n";
108 Exception exc(
"TGeoMaterialInterface::findNextBoundary ==> maximum number of iterations exceeded",__LINE__,__FILE__);
114 if (s + safety > fabs(sMax)) {
116 debugOut <<
" next boundary is further away than sMax \n";
117 return stepSign*(s + safety);
121 if (slDist < delta) {
123 debugOut <<
" very close to the boundary -> return"
124 <<
" stepSign*(s + slDist) = "
125 << stepSign <<
"*(" << s + slDist <<
")\n";
126 return stepSign*(s + slDist);
135 rep->
RKPropagate(state7,
nullptr, SA, stepSign*(s + step), varField);
139 double dist2 = (pow(state7[0] - oldState7[0], 2)
140 + pow(state7[1] - oldState7[1], 2)
141 + pow(state7[2] - oldState7[2], 2));
143 double maxDeviation2 = 0.25*(step*step - dist2);
146 && maxDeviation2 > epsilon*epsilon) {
155 step = std::max(step / 2, safety);
157 gGeoManager->PushPoint();
158 bool volChanged =
initTrack(state7[0], state7[1], state7[2],
159 stepSign*state7[3], stepSign*state7[4],
166 gGeoManager->PopPoint();
172 if (step <= safety) {
174 debugOut <<
" step <= safety, return stepSign*(s + step) = " << stepSign*(s +
step) <<
"\n";
175 return stepSign*(s +
step);
180 step = std::max(step / 2, safety);
186 gGeoManager->PopDummy();
188 gGeoManager->FindNextBoundary(fabs(sMax) - s);
189 safety = gGeoManager->GetSafeDistance();
190 slDist = gGeoManager->GetStep();
200 debugOut <<
" safety = " << safety <<
"; step, slDist = " << step <<
"\n";
216 1.e15, 19.2, 41.8, 40.0, 63.7, 76.0, 78., 82.0,
217 95.0, 115.0, 137.0, 149.0, 156.0, 166.0, 173.0, 173.0,
218 180.0, 174.0, 188.0, 190.0, 191.0, 216.0, 233.0, 245.0,
219 257.0, 272.0, 286.0, 297.0, 311.0, 322.0, 330.0, 334.0,
220 350.0, 347.0, 348.0, 343.0, 352.0, 363.0, 366.0, 379.0,
221 393.0, 417.0, 424.0, 428.0, 441.0, 449.0, 470.0, 470.0,
222 469.0, 488.0, 488.0, 487.0, 485.0, 491.0, 482.0, 488.0,
223 491.0, 501.0, 523.0, 535.0, 546.0, 560.0, 574.0, 580.0,
224 591.0, 614.0, 628.0, 650.0, 658.0, 674.0, 684.0, 694.0,
225 705.0, 718.0, 727.0, 736.0, 746.0, 757.0, 790.0, 790.0,
226 800.0, 810.0, 823.0, 823.0, 830.0, 825.0, 794.0, 827.0,
227 826.0, 841.0, 847.0, 878.0, 890.0, };
230 34.5388, 2.95491, 3.7329, 3.68888, 4.15418, 4.33073, 4.35671, 4.40672,
231 4.55388, 4.74493, 4.91998, 5.00395, 5.04986, 5.11199, 5.15329, 5.15329,
232 5.19296, 5.15906, 5.23644, 5.24702, 5.25227, 5.37528, 5.45104, 5.50126,
233 5.54908, 5.6058, 5.65599, 5.69373, 5.73979, 5.77455, 5.79909, 5.81114,
234 5.85793, 5.84932, 5.8522, 5.83773, 5.86363, 5.8944, 5.90263, 5.93754,
235 5.97381, 6.03309, 6.04973, 6.05912, 6.08904, 6.10702, 6.15273, 6.15273,
236 6.1506, 6.19032, 6.19032, 6.18826, 6.18415, 6.19644, 6.17794, 6.19032,
237 6.19644, 6.21661, 6.25958, 6.28227, 6.30262, 6.32794, 6.35263, 6.36303,
238 6.38182, 6.41999, 6.44254, 6.47697, 6.4892, 6.51323, 6.52796, 6.54247,
239 6.5582, 6.57647, 6.58893, 6.60123, 6.61473, 6.62936, 6.67203, 6.67203,
240 6.68461, 6.69703, 6.71296, 6.71296, 6.72143, 6.71538, 6.67708, 6.7178,
241 6.71659, 6.73459, 6.7417, 6.77765, 6.79122, };
256 if (mat->IsMixture()) {
261 TGeoMixture* mix = (TGeoMixture*)mat;
262 for (
int i = 0;
i < mix->GetNelements(); ++
i) {
263 int index = int(mix->GetZmixt()[
i]);
264 double weight = mix->GetWmixt()[
i] * mix->GetZmixt()[
i] / mix->GetAmixt()[
i];
273 int index = int(mat->GetZ());