GENFIT  Rev:NoNumberAvailable
TGeoMaterialInterface.cc
Go to the documentation of this file.
1 /* Copyright 2008-2014, Technische Universitaet Muenchen,
2  Authors: Christian Hoeppner & Sebastian Neubert & Johannes Rauch
3 
4  This file is part of GENFIT.
5 
6  GENFIT is free software: you can redistribute it and/or modify
7  it under the terms of the GNU Lesser General Public License as published
8  by the Free Software Foundation, either version 3 of the License, or
9  (at your option) any later version.
10 
11  GENFIT is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU Lesser General Public License for more details.
15 
16  You should have received a copy of the GNU Lesser General Public License
17  along with GENFIT. If not, see <http://www.gnu.org/licenses/>.
18 */
19 
20 #include "TGeoMaterialInterface.h"
21 #include "Exception.h"
22 #include "IO.h"
23 
24 #include <TGeoMedium.h>
25 #include <TGeoMaterial.h>
26 #include <TGeoManager.h>
27 #include <assert.h>
28 #include <math.h>
29 
30 
31 namespace genfit {
32 
33 double MeanExcEnergy_get(int Z);
34 double MeanExcEnergy_get(TGeoMaterial*);
35 
36 
37 bool
38 TGeoMaterialInterface::initTrack(double posX, double posY, double posZ,
39  double dirX, double dirY, double dirZ){
40  #ifdef DEBUG
41  debugOut << "TGeoMaterialInterface::initTrack. \n";
42  debugOut << "Pos "; TVector3(posX, posY, posZ).Print();
43  debugOut << "Dir "; TVector3(dirX, dirY, dirZ).Print();
44  #endif
45 
46  // Move to the new point.
47  bool result = !gGeoManager->IsSameLocation(posX, posY, posZ, kTRUE);
48  // Set the intended direction.
49  gGeoManager->SetCurrentDirection(dirX, dirY, dirZ);
50 
51  if (debugLvl_ > 0) {
52  debugOut << " TGeoMaterialInterface::initTrack at \n";
53  debugOut << " position: "; TVector3(posX, posY, posZ).Print();
54  debugOut << " direction: "; TVector3(dirX, dirY, dirZ).Print();
55  }
56 
57  return result;
58 }
59 
60 
61 void
63  double& Z,
64  double& A,
65  double& radiationLength,
66  double& mEE){
67 
68  TGeoMaterial* mat = gGeoManager->GetCurrentVolume()->GetMedium()->GetMaterial();
69 
70  density = mat->GetDensity();
71  Z = mat->GetZ();
72  A = mat->GetA();
73  radiationLength = mat->GetRadLen();
74  mEE = MeanExcEnergy_get(mat);
75 
76 }
77 
78 
79 void
81 
82  TGeoMaterial* mat = gGeoManager->GetCurrentVolume()->GetMedium()->GetMaterial();
83 
84  parameters.setMaterialProperties(mat->GetDensity(),
85  mat->GetZ(),
86  mat->GetA(),
87  mat->GetRadLen(),
88  MeanExcEnergy_get(mat));
89 
90 }
91 
92 
93 double
95  const M1x7& stateOrig,
96  double sMax, // signed
97  bool varField)
98 {
99  const double delta(1.E-2); // cm, distance limit beneath which straight-line steps are taken.
100  const double epsilon(1.E-1); // cm, allowed upper bound on arch
101  // deviation from straight line
102 
103  M1x3 SA;
104  M1x7 state7, oldState7;
105  oldState7 = stateOrig;
106 
107  int stepSign(sMax < 0 ? -1 : 1);
108 
109  double s = 0; // trajectory length to boundary
110 
111  const unsigned maxIt = 300;
112  unsigned it = 0;
113 
114  // Initialize the geometry to the current location (set by caller).
115  gGeoManager->FindNextBoundary(fabs(sMax) - s);
116  double safety = gGeoManager->GetSafeDistance(); // >= 0
117  double slDist = gGeoManager->GetStep();
118 
119  // this should not happen, but it happens sometimes.
120  // The reason are probably overlaps in the geometry.
121  // Without this "hack" many small steps with size of minStep will be made,
122  // until eventually the maximum number of iterations are exceeded and the extrapolation fails
123  if (slDist < safety)
124  slDist = safety;
125  double step = slDist;
126 
127  if (debugLvl_ > 0)
128  debugOut << " safety = " << safety << "; step, slDist = " << step << "\n";
129 
130  while (1) {
131  if (++it > maxIt){
132  Exception exc("TGeoMaterialInterface::findNextBoundary ==> maximum number of iterations exceeded",__LINE__,__FILE__);
133  exc.setFatal();
134  throw exc;
135  }
136 
137  // No boundary in sight?
138  if (s + safety > fabs(sMax)) {
139  if (debugLvl_ > 0)
140  debugOut << " next boundary is further away than sMax \n";
141  return stepSign*(s + safety); //sMax;
142  }
143 
144  // Are we at the boundary?
145  if (slDist < delta) {
146  if (debugLvl_ > 0)
147  debugOut << " very close to the boundary -> return"
148  << " stepSign*(s + slDist) = "
149  << stepSign << "*(" << s + slDist << ")\n";
150  return stepSign*(s + slDist);
151  }
152 
153  // We have to find whether there's any boundary on our path.
154 
155  // Follow curved arch, then see if we may have missed a boundary.
156  // Always propagate complete way from original start to avoid
157  // inconsistent extrapolations.
158  state7 = stateOrig;
159  rep->RKPropagate(state7, NULL, SA, stepSign*(s + step), varField);
160 
161  // Straight line distance² between extrapolation finish and
162  // the end of the previously determined safe segment.
163  double dist2 = (pow(state7[0] - oldState7[0], 2)
164  + pow(state7[1] - oldState7[1], 2)
165  + pow(state7[2] - oldState7[2], 2));
166  // Maximal lateral deviation².
167  double maxDeviation2 = 0.25*(step*step - dist2);
168 
169  if (step > safety
170  && maxDeviation2 > epsilon*epsilon) {
171  // Need to take a shorter step to reliably estimate material,
172  // but only if we didn't move by safety. In order to avoid
173  // issues with roundoff we check 'step > safety' instead of
174  // 'step != safety'. If we moved by safety, there couldn't have
175  // been a boundary that we missed along the path, but we could
176  // be on the boundary.
177 
178  // Take a shorter step, but never shorter than safety.
179  step = std::max(step / 2, safety);
180  } else {
181  gGeoManager->PushPoint();
182  bool volChanged = initTrack(state7[0], state7[1], state7[2],
183  stepSign*state7[3], stepSign*state7[4],
184  stepSign*state7[5]);
185 
186  if (volChanged) {
187  if (debugLvl_ > 0)
188  debugOut << " volChanged\n";
189  // Move back to start.
190  gGeoManager->PopPoint();
191 
192  // Extrapolation may not take the exact step length we asked
193  // for, so it can happen that a requested step < safety takes
194  // us across the boundary. This is then the best estimate we
195  // can get of the distance to the boundary with the stepper.
196  if (step <= safety) {
197  if (debugLvl_ > 0)
198  debugOut << " step <= safety, return stepSign*(s + step) = " << stepSign*(s + step) << "\n";
199  return stepSign*(s + step);
200  }
201 
202  // Volume changed during the extrapolation. Take a shorter
203  // step, but never shorter than safety.
204  step = std::max(step / 2, safety);
205  } else {
206  // we're in the new place, the step was safe, advance
207  s += step;
208 
209  oldState7 = state7;
210  gGeoManager->PopDummy(); // Pop stack, but stay in place.
211 
212  gGeoManager->FindNextBoundary(fabs(sMax) - s);
213  safety = gGeoManager->GetSafeDistance();
214  slDist = gGeoManager->GetStep();
215 
216  // this should not happen, but it happens sometimes.
217  // The reason are probably overlaps in the geometry.
218  // Without this "hack" many small steps with size of minStep will be made,
219  // until eventually the maximum number of iterations are exceeded and the extrapolation fails
220  if (slDist < safety)
221  slDist = safety;
222  step = slDist;
223  if (debugLvl_ > 0)
224  debugOut << " safety = " << safety << "; step, slDist = " << step << "\n";
225  }
226  }
227  }
228 }
229 
230 
231 /*
232 Reference for elemental mean excitation energies at:
233 http://physics.nist.gov/PhysRefData/XrayMassCoef/tab1.html
234 
235 Code ported from GEANT 3
236 */
237 
238 const int MeanExcEnergy_NELEMENTS = 93; // 0 = vacuum, 1 = hydrogen, 92 = uranium
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, };
252 // Logarithms of the previous table, used to calculate mixtures.
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, };
266 
267 
268 double
269 MeanExcEnergy_get(int Z, bool logs = false) {
270  assert(Z >= 0 && Z < MeanExcEnergy_NELEMENTS);
271  if (logs)
272  return MeanExcEnergy_logs[Z];
273  else
274  return MeanExcEnergy_vals[Z];
275 }
276 
277 
278 double
279 MeanExcEnergy_get(TGeoMaterial* mat) {
280  if (mat->IsMixture()) {
281  // The mean excitation energy is calculated as the geometric mean
282  // of the mEEs of the components, weighted for abundance.
283  double logMEE = 0.;
284  double denom = 0.;
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];
289  logMEE += weight * MeanExcEnergy_get(index, true);
290  denom += weight;
291  }
292  logMEE /= denom;
293  return exp(logMEE);
294  }
295 
296  // not a mixture
297  int index = int(mat->GetZ());
298  return MeanExcEnergy_get(index, false);
299 }
300 
301 
302 } /* End of namespace genfit */
genfit::TGeoMaterialInterface::getMaterialParameters
void getMaterialParameters(double &density, double &Z, double &A, double &radiationLength, double &mEE)
Get material parameters in current material.
Definition: TGeoMaterialInterface.cc:62
genfit::AbsMaterialInterface::debugLvl_
unsigned int debugLvl_
Definition: AbsMaterialInterface.h:75
genfit::RKTrackRep::RKPropagate
double RKPropagate(M1x7 &state7, M7x7 *jacobian, M1x3 &SA, double S, bool varField=true, bool calcOnlyLastRowOfJ=false) const
The actual Runge Kutta propagation.
Definition: RKTrackRep.cc:1272
genfit::RKMatrix< 1, 7 >
TGeoMaterialInterface.h
genfit::MeanExcEnergy_NELEMENTS
const int MeanExcEnergy_NELEMENTS
Definition: TGeoMaterialInterface.cc:238
genfit::MaterialProperties::setMaterialProperties
void setMaterialProperties(const double &density, const double &Z, const double &A, const double &radiationLength, const double &mEE)
Definition: MaterialProperties.cc:57
genfit::TGeoMaterialInterface::findNextBoundary
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....
Definition: TGeoMaterialInterface.cc:94
genfit::Exception
Exception class for error handling in GENFIT (provides storage for diagnostic information)
Definition: Exception.h:48
genfit
Defines for I/O streams used for error and debug printing.
Definition: AbsFinitePlane.cc:22
genfit::RKTrackRep
AbsTrackRep with 5D track parameterization in plane coordinates: (q/p, u', v', u, v)
Definition: RKTrackRep.h:71
genfit::MeanExcEnergy_get
double MeanExcEnergy_get(int Z)
genfit::Exception::setFatal
void setFatal(bool b=true)
Set fatal flag.
Definition: Exception.h:61
IO.h
genfit::MeanExcEnergy_logs
const double MeanExcEnergy_logs[MeanExcEnergy_NELEMENTS]
Definition: TGeoMaterialInterface.cc:253
genfit::MeanExcEnergy_vals
const double MeanExcEnergy_vals[MeanExcEnergy_NELEMENTS]
Definition: TGeoMaterialInterface.cc:239
Exception.h
genfit::TGeoMaterialInterface::initTrack
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...
Definition: TGeoMaterialInterface.cc:38
genfit::MaterialProperties
Material properties needed e.g. for material effects calculation.
Definition: MaterialProperties.h:35
genfit::debugOut
std::ostream debugOut