Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TGeoMaterialInterface.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TGeoMaterialInterface.cc
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 
62 
63  TGeoMaterial* mat = gGeoManager->GetCurrentVolume()->GetMedium()->GetMaterial();
64  return Material(mat->GetDensity(), mat->GetZ(), mat->GetA(), mat->GetRadLen(), MeanExcEnergy_get(mat));
65 
66 }
67 
68 
69 double
71  const M1x7& stateOrig,
72  double sMax, // signed
73  bool varField)
74 {
75  const double delta(1.E-2); // cm, distance limit beneath which straight-line steps are taken.
76  const double epsilon(1.E-1); // cm, allowed upper bound on arch
77  // deviation from straight line
78 
79  M1x3 SA;
80  M1x7 state7, oldState7;
81  oldState7 = stateOrig;
82 
83  int stepSign(sMax < 0 ? -1 : 1);
84 
85  double s = 0; // trajectory length to boundary
86 
87  const unsigned maxIt = 300;
88  unsigned it = 0;
89 
90  // Initialize the geometry to the current location (set by caller).
91  gGeoManager->FindNextBoundary(fabs(sMax) - s);
92  double safety = gGeoManager->GetSafeDistance(); // >= 0
93  double slDist = gGeoManager->GetStep();
94 
95  // this should not happen, but it happens sometimes.
96  // The reason are probably overlaps in the geometry.
97  // Without this "hack" many small steps with size of minStep will be made,
98  // until eventually the maximum number of iterations are exceeded and the extrapolation fails
99  if (slDist < safety)
100  slDist = safety;
101  double step = slDist;
102 
103  if (debugLvl_ > 0)
104  debugOut << " safety = " << safety << "; step, slDist = " << step << "\n";
105 
106  while (1) {
107  if (++it > maxIt){
108  Exception exc("TGeoMaterialInterface::findNextBoundary ==> maximum number of iterations exceeded",__LINE__,__FILE__);
109  exc.setFatal();
110  throw exc;
111  }
112 
113  // No boundary in sight?
114  if (s + safety > fabs(sMax)) {
115  if (debugLvl_ > 0)
116  debugOut << " next boundary is further away than sMax \n";
117  return stepSign*(s + safety); //sMax;
118  }
119 
120  // Are we at the boundary?
121  if (slDist < delta) {
122  if (debugLvl_ > 0)
123  debugOut << " very close to the boundary -> return"
124  << " stepSign*(s + slDist) = "
125  << stepSign << "*(" << s + slDist << ")\n";
126  return stepSign*(s + slDist);
127  }
128 
129  // We have to find whether there's any boundary on our path.
130 
131  // Follow curved arch, then see if we may have missed a boundary.
132  // Always propagate complete way from original start to avoid
133  // inconsistent extrapolations.
134  state7 = stateOrig;
135  rep->RKPropagate(state7, nullptr, SA, stepSign*(s + step), varField);
136 
137  // Straight line distance² between extrapolation finish and
138  // the end of the previously determined safe segment.
139  double dist2 = (pow(state7[0] - oldState7[0], 2)
140  + pow(state7[1] - oldState7[1], 2)
141  + pow(state7[2] - oldState7[2], 2));
142  // Maximal lateral deviation².
143  double maxDeviation2 = 0.25*(step*step - dist2);
144 
145  if (step > safety
146  && maxDeviation2 > epsilon*epsilon) {
147  // Need to take a shorter step to reliably estimate material,
148  // but only if we didn't move by safety. In order to avoid
149  // issues with roundoff we check 'step > safety' instead of
150  // 'step != safety'. If we moved by safety, there couldn't have
151  // been a boundary that we missed along the path, but we could
152  // be on the boundary.
153 
154  // Take a shorter step, but never shorter than safety.
155  step = std::max(step / 2, safety);
156  } else {
157  gGeoManager->PushPoint();
158  bool volChanged = initTrack(state7[0], state7[1], state7[2],
159  stepSign*state7[3], stepSign*state7[4],
160  stepSign*state7[5]);
161 
162  if (volChanged) {
163  if (debugLvl_ > 0)
164  debugOut << " volChanged\n";
165  // Move back to start.
166  gGeoManager->PopPoint();
167 
168  // Extrapolation may not take the exact step length we asked
169  // for, so it can happen that a requested step < safety takes
170  // us across the boundary. This is then the best estimate we
171  // can get of the distance to the boundary with the stepper.
172  if (step <= safety) {
173  if (debugLvl_ > 0)
174  debugOut << " step <= safety, return stepSign*(s + step) = " << stepSign*(s + step) << "\n";
175  return stepSign*(s + step);
176  }
177 
178  // Volume changed during the extrapolation. Take a shorter
179  // step, but never shorter than safety.
180  step = std::max(step / 2, safety);
181  } else {
182  // we're in the new place, the step was safe, advance
183  s += step;
184 
185  oldState7 = state7;
186  gGeoManager->PopDummy(); // Pop stack, but stay in place.
187 
188  gGeoManager->FindNextBoundary(fabs(sMax) - s);
189  safety = gGeoManager->GetSafeDistance();
190  slDist = gGeoManager->GetStep();
191 
192  // this should not happen, but it happens sometimes.
193  // The reason are probably overlaps in the geometry.
194  // Without this "hack" many small steps with size of minStep will be made,
195  // until eventually the maximum number of iterations are exceeded and the extrapolation fails
196  if (slDist < safety)
197  slDist = safety;
198  step = slDist;
199  if (debugLvl_ > 0)
200  debugOut << " safety = " << safety << "; step, slDist = " << step << "\n";
201  }
202  }
203  }
204 }
205 
206 
207 /*
208 Reference for elemental mean excitation energies at:
209 http://physics.nist.gov/PhysRefData/XrayMassCoef/tab1.html
210 
211 Code ported from GEANT 3
212 */
213 
214 const int MeanExcEnergy_NELEMENTS = 93; // 0 = vacuum, 1 = hydrogen, 92 = uranium
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, };
228 // Logarithms of the previous table, used to calculate mixtures.
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, };
242 
243 
244 double
245 MeanExcEnergy_get(int Z, bool logs = false) {
246  assert(Z >= 0 && Z < MeanExcEnergy_NELEMENTS);
247  if (logs)
248  return MeanExcEnergy_logs[Z];
249  else
250  return MeanExcEnergy_vals[Z];
251 }
252 
253 
254 double
255 MeanExcEnergy_get(TGeoMaterial* mat) {
256  if (mat->IsMixture()) {
257  // The mean excitation energy is calculated as the geometric mean
258  // of the mEEs of the components, weighted for abundance.
259  double logMEE = 0.;
260  double denom = 0.;
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];
265  logMEE += weight * MeanExcEnergy_get(index, true);
266  denom += weight;
267  }
268  logMEE /= denom;
269  return exp(logMEE);
270  }
271 
272  // not a mixture
273  int index = int(mat->GetZ());
274  return MeanExcEnergy_get(index, false);
275 }
276 
277 
278 } /* End of namespace genfit */