trackReps/src/TGeoMaterialInterface.cc

0001 /* Copyright 2008-2014, Technische Universitaet Muenchen,
0002    Authors: Christian Hoeppner & Sebastian Neubert & Johannes Rauch
0003
0004    This file is part of GENFIT.
0005
0006    GENFIT is free software: you can redistribute it and/or modify
0007    it under the terms of the GNU Lesser General Public License as published
0008    by the Free Software Foundation, either version 3 of the License, or
0009    (at your option) any later version.
0010
0011    GENFIT is distributed in the hope that it will be useful,
0012    but WITHOUT ANY WARRANTY; without even the implied warranty of
0013    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0014    GNU Lesser General Public License for more details.
0015
0016    You should have received a copy of the GNU Lesser General Public License
0017    along with GENFIT.  If not, see <http://www.gnu.org/licenses/>.
0018 */
0019
0020 #include "TGeoMaterialInterface.h"
0021 #include "Exception.h"
0022 #include "IO.h"
0023
0024 #include <TGeoMedium.h>
0025 #include <TGeoMaterial.h>
0026 #include <TGeoManager.h>
0027 #include <assert.h>
0028 #include <math.h>
0029
0030
0031 namespace genfit {
0032
0033 double MeanExcEnergy_get(int Z);
0034 double MeanExcEnergy_get(TGeoMaterial*);
0035
0036
0037 bool
0038 TGeoMaterialInterface::initTrack(double posX, double posY, double posZ,
0039                                    double dirX, double dirY, double dirZ){
0040   #ifdef DEBUG
0041   debugOut << "TGeoMaterialInterface::initTrack. \n";
0042   debugOut << "Pos    "; TVector3(posX, posY, posZ).Print();
0043   debugOut << "Dir    "; TVector3(dirX, dirY, dirZ).Print();
0044   #endif
0045
0046   // Move to the new point.
0047   bool result = !gGeoManager->IsSameLocation(posX, posY, posZ, kTRUE);
0048   // Set the intended direction.
0049   gGeoManager->SetCurrentDirection(dirX, dirY, dirZ);
0050
0051   if (debugLvl_ > 0) {
0052     debugOut << "      TGeoMaterialInterface::initTrack at \n";
0053     debugOut << "      position:  "; TVector3(posX, posY, posZ).Print();
0054     debugOut << "      direction: "; TVector3(dirX, dirY, dirZ).Print();
0055   }
0056
0057   return result;
0058 }
0059
0060
0061 Material TGeoMaterialInterface::getMaterialParameters() {
0062
0063   TGeoMaterial* mat = gGeoManager->GetCurrentVolume()->GetMedium()->GetMaterial();
0064   return Material(mat->GetDensity(), mat->GetZ(), mat->GetA(), mat->GetRadLen(), MeanExcEnergy_get(mat));
0065
0066 }
0067
0068
0069 double
0070 TGeoMaterialInterface::findNextBoundary(const RKTrackRep* rep,
0071                                           const M1x7& stateOrig,
0072                                           double sMax, // signed
0073                                           bool varField)
0074 {
0075   const double delta(1.E-2); // cm, distance limit beneath which straight-line steps are taken.
0076   const double epsilon(1.E-1); // cm, allowed upper bound on arch
0077   // deviation from straight line
0078
0079   M1x3 SA;
0080   M1x7 state7, oldState7;
0081   oldState7 = stateOrig;
0082
0083   int stepSign(sMax < 0 ? -1 : 1);
0084
0085   double s = 0;  // trajectory length to boundary
0086
0087   const unsigned maxIt = 300;
0088   unsigned it = 0;
0089
0090   // Initialize the geometry to the current location (set by caller).
0091   gGeoManager->FindNextBoundary(fabs(sMax) - s);
0092   double safety = gGeoManager->GetSafeDistance(); // >= 0
0093   double slDist = gGeoManager->GetStep();
0094
0095   // this should not happen, but it happens sometimes.
0096   // The reason are probably overlaps in the geometry.
0097   // Without this "hack" many small steps with size of minStep will be made,
0098   // until eventually the maximum number of iterations are exceeded and the extrapolation fails
0099   if (slDist < safety)
0100     slDist = safety;
0101   double step = slDist;
0102
0103   if (debugLvl_ > 0)
0104     debugOut << "   safety = " << safety << "; step, slDist = " << step << "\n";
0105
0106   while (1) {
0107     if (++it > maxIt){
0108       Exception exc("TGeoMaterialInterface::findNextBoundary ==> maximum number of iterations exceeded",__LINE__,__FILE__);
0109       exc.setFatal();
0110       throw exc;
0111     }
0112
0113     // No boundary in sight?
0114     if (s + safety > fabs(sMax)) {
0115       if (debugLvl_ > 0)
0116         debugOut << "   next boundary is further away than sMax \n";
0117       return stepSign*(s + safety); //sMax;
0118     }
0119
0120     // Are we at the boundary?
0121     if (slDist < delta) {
0122       if (debugLvl_ > 0)
0123         debugOut << "   very close to the boundary -> return"
0124         << " stepSign*(s + slDist) = "
0125         << stepSign << "*(" << s + slDist << ")\n";
0126       return stepSign*(s + slDist);
0127     }
0128
0129     // We have to find whether there's any boundary on our path.
0130
0131     // Follow curved arch, then see if we may have missed a boundary.
0132     // Always propagate complete way from original start to avoid
0133     // inconsistent extrapolations.
0134     state7 = stateOrig;
0135     rep->RKPropagate(state7, nullptr, SA, stepSign*(s + step), varField);
0136
0137     // Straight line distance² between extrapolation finish and
0138     // the end of the previously determined safe segment.
0139     double dist2 = (pow(state7[0] - oldState7[0], 2)
0140         + pow(state7[1] - oldState7[1], 2)
0141         + pow(state7[2] - oldState7[2], 2));
0142     // Maximal lateral deviation².
0143     double maxDeviation2 = 0.25*(step*step - dist2);
0144
0145     if (step > safety
0146         && maxDeviation2 > epsilon*epsilon) {
0147       // Need to take a shorter step to reliably estimate material,
0148       // but only if we didn't move by safety.  In order to avoid
0149       // issues with roundoff we check 'step > safety' instead of
0150       // 'step != safety'.  If we moved by safety, there couldn't have
0151       // been a boundary that we missed along the path, but we could
0152       // be on the boundary.
0153
0154       // Take a shorter step, but never shorter than safety.
0155       step = std::max(step / 2, safety);
0156     } else {
0157       gGeoManager->PushPoint();
0158       bool volChanged = initTrack(state7[0], state7[1], state7[2],
0159           stepSign*state7[3], stepSign*state7[4],
0160           stepSign*state7[5]);
0161
0162       if (volChanged) {
0163         if (debugLvl_ > 0)
0164           debugOut << "   volChanged\n";
0165         // Move back to start.
0166         gGeoManager->PopPoint();
0167
0168         // Extrapolation may not take the exact step length we asked
0169         // for, so it can happen that a requested step < safety takes
0170         // us across the boundary.  This is then the best estimate we
0171         // can get of the distance to the boundary with the stepper.
0172         if (step <= safety) {
0173           if (debugLvl_ > 0)
0174             debugOut << "   step <= safety, return stepSign*(s + step) = " << stepSign*(s + step) << "\n";
0175           return stepSign*(s + step);
0176         }
0177
0178         // Volume changed during the extrapolation.  Take a shorter
0179         // step, but never shorter than safety.
0180         step = std::max(step / 2, safety);
0181       } else {
0182         // we're in the new place, the step was safe, advance
0183         s += step;
0184
0185         oldState7 = state7;
0186         gGeoManager->PopDummy();  // Pop stack, but stay in place.
0187
0188         gGeoManager->FindNextBoundary(fabs(sMax) - s);
0189         safety = gGeoManager->GetSafeDistance();
0190         slDist = gGeoManager->GetStep();
0191
0192         // this should not happen, but it happens sometimes.
0193         // The reason are probably overlaps in the geometry.
0194         // Without this "hack" many small steps with size of minStep will be made,
0195         // until eventually the maximum number of iterations are exceeded and the extrapolation fails
0196         if (slDist < safety)
0197           slDist = safety;
0198         step = slDist;
0199         if (debugLvl_ > 0)
0200           debugOut << "   safety = " << safety << "; step, slDist = " << step << "\n";
0201       }
0202     }
0203   }
0204 }
0205
0206
0207 /*
0208 Reference for elemental mean excitation energies at:
0209 http://physics.nist.gov/PhysRefData/XrayMassCoef/tab1.html
0210
0211 Code ported from GEANT 3
0212 */
0213
0214 const int MeanExcEnergy_NELEMENTS = 93; // 0 = vacuum, 1 = hydrogen, 92 = uranium
0215 const double MeanExcEnergy_vals[MeanExcEnergy_NELEMENTS] = {
0216   1.e15, 19.2, 41.8, 40.0, 63.7, 76.0, 78., 82.0,
0217   95.0, 115.0, 137.0, 149.0, 156.0, 166.0, 173.0, 173.0,
0218   180.0, 174.0, 188.0, 190.0, 191.0, 216.0, 233.0, 245.0,
0219   257.0, 272.0, 286.0, 297.0, 311.0, 322.0, 330.0, 334.0,
0220   350.0, 347.0, 348.0, 343.0, 352.0, 363.0, 366.0, 379.0,
0221   393.0, 417.0, 424.0, 428.0, 441.0, 449.0, 470.0, 470.0,
0222   469.0, 488.0, 488.0, 487.0, 485.0, 491.0, 482.0, 488.0,
0223   491.0, 501.0, 523.0, 535.0, 546.0, 560.0, 574.0, 580.0,
0224   591.0, 614.0, 628.0, 650.0, 658.0, 674.0, 684.0, 694.0,
0225   705.0, 718.0, 727.0, 736.0, 746.0, 757.0, 790.0, 790.0,
0226   800.0, 810.0, 823.0, 823.0, 830.0, 825.0, 794.0, 827.0,
0227   826.0, 841.0, 847.0, 878.0, 890.0, };
0228 // Logarithms of the previous table, used to calculate mixtures.
0229 const double MeanExcEnergy_logs[MeanExcEnergy_NELEMENTS] = {
0230   34.5388, 2.95491, 3.7329, 3.68888, 4.15418, 4.33073, 4.35671, 4.40672,
0231   4.55388, 4.74493, 4.91998, 5.00395, 5.04986, 5.11199, 5.15329, 5.15329,
0232   5.19296, 5.15906, 5.23644, 5.24702, 5.25227, 5.37528, 5.45104, 5.50126,
0233   5.54908, 5.6058, 5.65599, 5.69373, 5.73979, 5.77455, 5.79909, 5.81114,
0234   5.85793, 5.84932, 5.8522, 5.83773, 5.86363, 5.8944, 5.90263, 5.93754,
0235   5.97381, 6.03309, 6.04973, 6.05912, 6.08904, 6.10702, 6.15273, 6.15273,
0236   6.1506, 6.19032, 6.19032, 6.18826, 6.18415, 6.19644, 6.17794, 6.19032,
0237   6.19644, 6.21661, 6.25958, 6.28227, 6.30262, 6.32794, 6.35263, 6.36303,
0238   6.38182, 6.41999, 6.44254, 6.47697, 6.4892, 6.51323, 6.52796, 6.54247,
0239   6.5582, 6.57647, 6.58893, 6.60123, 6.61473, 6.62936, 6.67203, 6.67203,
0240   6.68461, 6.69703, 6.71296, 6.71296, 6.72143, 6.71538, 6.67708, 6.7178,
0241   6.71659, 6.73459, 6.7417, 6.77765, 6.79122, };
0242
0243
0244 double
0245 MeanExcEnergy_get(int Z, bool logs = false) {
0246   assert(Z >= 0 && Z < MeanExcEnergy_NELEMENTS);
0247   if (logs)
0248     return MeanExcEnergy_logs[Z];
0249   else
0250     return MeanExcEnergy_vals[Z];
0251 }
0252
0253
0254 double
0255 MeanExcEnergy_get(TGeoMaterial* mat) {
0256   if (mat->IsMixture()) {
0257     // The mean excitation energy is calculated as the geometric mean
0258     // of the mEEs of the components, weighted for abundance.
0259     double logMEE = 0.;
0260     double denom  = 0.;
0261     TGeoMixture* mix = (TGeoMixture*)mat;
0262     for (int i = 0; i < mix->GetNelements(); ++i) {
0263       int index = int(mix->GetZmixt()[i]);
0264       double weight = mix->GetWmixt()[i] * mix->GetZmixt()[i] / mix->GetAmixt()[i];
0265       logMEE += weight * MeanExcEnergy_get(index, true);
0266       denom  += weight;
0267     }
0268     logMEE /= denom;
0269     return exp(logMEE);
0270   }
0271
0272   // not a mixture
0273   int index = int(mat->GetZ());
0274   return MeanExcEnergy_get(index, false);
0275 }
0276
0277
0278 } /* End of namespace genfit */