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();
65 double& radiationLength,
68 TGeoMaterial* mat = gGeoManager->GetCurrentVolume()->GetMedium()->GetMaterial();
70 density = mat->GetDensity();
73 radiationLength = mat->GetRadLen();
82 TGeoMaterial* mat = gGeoManager->GetCurrentVolume()->GetMedium()->GetMaterial();
95 const M1x7& stateOrig,
99 const double delta(1.E-2);
100 const double epsilon(1.E-1);
104 M1x7 state7, oldState7;
105 oldState7 = stateOrig;
107 int stepSign(sMax < 0 ? -1 : 1);
111 const unsigned maxIt = 300;
115 gGeoManager->FindNextBoundary(fabs(sMax) - s);
116 double safety = gGeoManager->GetSafeDistance();
117 double slDist = gGeoManager->GetStep();
125 double step = slDist;
128 debugOut <<
" safety = " << safety <<
"; step, slDist = " << step <<
"\n";
132 Exception exc(
"TGeoMaterialInterface::findNextBoundary ==> maximum number of iterations exceeded",__LINE__,__FILE__);
138 if (s + safety > fabs(sMax)) {
140 debugOut <<
" next boundary is further away than sMax \n";
141 return stepSign*(s + safety);
145 if (slDist < delta) {
147 debugOut <<
" very close to the boundary -> return"
148 <<
" stepSign*(s + slDist) = "
149 << stepSign <<
"*(" << s + slDist <<
")\n";
150 return stepSign*(s + slDist);
159 rep->
RKPropagate(state7, NULL, SA, stepSign*(s + step), varField);
163 double dist2 = (pow(state7[0] - oldState7[0], 2)
164 + pow(state7[1] - oldState7[1], 2)
165 + pow(state7[2] - oldState7[2], 2));
167 double maxDeviation2 = 0.25*(step*step - dist2);
170 && maxDeviation2 > epsilon*epsilon) {
179 step = std::max(step / 2, safety);
181 gGeoManager->PushPoint();
182 bool volChanged =
initTrack(state7[0], state7[1], state7[2],
183 stepSign*state7[3], stepSign*state7[4],
190 gGeoManager->PopPoint();
196 if (step <= safety) {
198 debugOut <<
" step <= safety, return stepSign*(s + step) = " << stepSign*(s + step) <<
"\n";
199 return stepSign*(s + step);
204 step = std::max(step / 2, safety);
210 gGeoManager->PopDummy();
212 gGeoManager->FindNextBoundary(fabs(sMax) - s);
213 safety = gGeoManager->GetSafeDistance();
214 slDist = gGeoManager->GetStep();
224 debugOut <<
" safety = " << safety <<
"; step, slDist = " << step <<
"\n";
240 1.e15, 19.2, 41.8, 40.0, 63.7, 76.0, 78., 82.0,
241 95.0, 115.0, 137.0, 149.0, 156.0, 166.0, 173.0, 173.0,
242 180.0, 174.0, 188.0, 190.0, 191.0, 216.0, 233.0, 245.0,
243 257.0, 272.0, 286.0, 297.0, 311.0, 322.0, 330.0, 334.0,
244 350.0, 347.0, 348.0, 343.0, 352.0, 363.0, 366.0, 379.0,
245 393.0, 417.0, 424.0, 428.0, 441.0, 449.0, 470.0, 470.0,
246 469.0, 488.0, 488.0, 487.0, 485.0, 491.0, 482.0, 488.0,
247 491.0, 501.0, 523.0, 535.0, 546.0, 560.0, 574.0, 580.0,
248 591.0, 614.0, 628.0, 650.0, 658.0, 674.0, 684.0, 694.0,
249 705.0, 718.0, 727.0, 736.0, 746.0, 757.0, 790.0, 790.0,
250 800.0, 810.0, 823.0, 823.0, 830.0, 825.0, 794.0, 827.0,
251 826.0, 841.0, 847.0, 878.0, 890.0, };
254 34.5388, 2.95491, 3.7329, 3.68888, 4.15418, 4.33073, 4.35671, 4.40672,
255 4.55388, 4.74493, 4.91998, 5.00395, 5.04986, 5.11199, 5.15329, 5.15329,
256 5.19296, 5.15906, 5.23644, 5.24702, 5.25227, 5.37528, 5.45104, 5.50126,
257 5.54908, 5.6058, 5.65599, 5.69373, 5.73979, 5.77455, 5.79909, 5.81114,
258 5.85793, 5.84932, 5.8522, 5.83773, 5.86363, 5.8944, 5.90263, 5.93754,
259 5.97381, 6.03309, 6.04973, 6.05912, 6.08904, 6.10702, 6.15273, 6.15273,
260 6.1506, 6.19032, 6.19032, 6.18826, 6.18415, 6.19644, 6.17794, 6.19032,
261 6.19644, 6.21661, 6.25958, 6.28227, 6.30262, 6.32794, 6.35263, 6.36303,
262 6.38182, 6.41999, 6.44254, 6.47697, 6.4892, 6.51323, 6.52796, 6.54247,
263 6.5582, 6.57647, 6.58893, 6.60123, 6.61473, 6.62936, 6.67203, 6.67203,
264 6.68461, 6.69703, 6.71296, 6.71296, 6.72143, 6.71538, 6.67708, 6.7178,
265 6.71659, 6.73459, 6.7417, 6.77765, 6.79122, };
280 if (mat->IsMixture()) {
285 TGeoMixture* mix = (TGeoMixture*)mat;
286 for (
int i = 0; i < mix->GetNelements(); ++i) {
287 int index = int(mix->GetZmixt()[i]);
288 double weight = mix->GetWmixt()[i] * mix->GetZmixt()[i] / mix->GetAmixt()[i];
297 int index = int(mat->GetZ());