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, };
270 assert(Z >= 0 && Z < MeanExcEnergy_NELEMENTS);
272 return MeanExcEnergy_logs[Z];
274 return MeanExcEnergy_vals[Z];
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());
double RKPropagate(M1x7 &state7, M7x7 *jacobian, M1x3 &SA, double S, bool varField=true, bool calcOnlyLastRowOfJ=false) const
The actual Runge Kutta propagation.
double findNextBoundary(const RKTrackRep *rep, const M1x7 &state7, double sMax, bool varField=true)
Make a step (following the curvature) until step length sMax or the next boundary is reached...
const int MeanExcEnergy_NELEMENTS
const double MeanExcEnergy_vals[MeanExcEnergy_NELEMENTS]
void setFatal(bool b=true)
Set fatal flag.
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
const double MeanExcEnergy_logs[MeanExcEnergy_NELEMENTS]
AbsTrackRep with 5D track parameterization in plane coordinates: (q/p, u', v', u, v) ...
Material properties needed e.g. for material effects calculation.
void setMaterialProperties(const double &density, const double &Z, const double &A, const double &radiationLength, const double &mEE)
double MeanExcEnergy_get(int Z)
bool initTrack(double posX, double posY, double posZ, double dirX, double dirY, double dirZ)
Initialize the navigator at given position and with given direction. Returns true if the volume chang...
void getMaterialParameters(double &density, double &Z, double &A, double &radiationLength, double &mEE)
Get material parameters in current material.
Defines for I/O streams used for error and debug printing.