include/G4HepEm/G4HepEmRunUtils.icc

0001
0002 #include "G4HepEmRunUtils.hh"
0003
0004 #include "G4HepEmMath.hh"
0005
0006 #include <cmath>
0007 #include <algorithm>
0008
0009
0010 // Roate the direction [u,v,w] given in the scattering frame to the lab frame.
0011 // Details: scattering is described relative to the [0,0,1] direction (i.e. scattering
0012 // frame). Therefore, after the new direction is computed relative to this [0,0,1]
0013 // original direction, the real original direction [u1,u2,u3] in the lab frame
0014 // needs to be accounted and the final new direction, i.e. in the lab frame is
0015 // computed.
0016 void RotateToReferenceFrame(double &u, double &v, double &w, const double* refDir) {
0017   double up = refDir[0]*refDir[0] + refDir[1]*refDir[1];
0018   if (up>0.) {
0019     up = std::sqrt(up);
0020     const double px = u;
0021     const double py = v;
0022     const double pz = w;
0023     u = (refDir[0]*refDir[2]*px - refDir[1]*py)/up + refDir[0]*pz;
0024     v = (refDir[1]*refDir[2]*px + refDir[0]*py)/up + refDir[1]*pz;
0025     w =    -up*px +             refDir[2]*pz;
0026   } else if (refDir[2]<0.) {       // phi=0  teta=pi
0027     u = -u;
0028     w = -w;
0029   }
0030 }
0031
0032 void RotateToReferenceFrame(double* dir, const double* refDir) {
0033   double up = refDir[0]*refDir[0] + refDir[1]*refDir[1];
0034   if (up>0.) {
0035     up = std::sqrt(up);
0036     const double px = dir[0];
0037     const double py = dir[1];
0038     const double pz = dir[2];
0039     dir[0] = (refDir[0]*refDir[2]*px - refDir[1]*py)/up + refDir[0]*pz;
0040     dir[1] = (refDir[1]*refDir[2]*px + refDir[0]*py)/up + refDir[1]*pz;
0041     dir[2] =    -up*px +             refDir[2]*pz;
0042   } else if (refDir[2]<0.) {       // phi=0  teta=pi
0043     dir[0] = -dir[0];
0044     dir[2] = -dir[2];
0045   }
0046 }
0047
0048 // use the improved, robust spline interpolation that I put in G4 10.6
0049 double GetSpline(double x1, double x2, double y1, double y2, double secderiv1, double secderiv2, double x)
0050 {
0051   // Unchecked precondition: x1 < x < x2
0052   const double dl = x2 - x1;
0053   // note: all corner cases of the previous methods are covered and eventually
0054   //       gives b=0/1 that results in y=y0\y_{N-1} if e<=x[0]/e>=x[N-1] or
0055   //       y=y_i/y_{i+1} if e<x[i]/e>=x[i+1] due to small numerical errors
0056   const double  b = G4HepEmMax(0., G4HepEmMin(1., (x - x1)/dl));
0057   const double os = 0.166666666667; // 1./6.
0058   const double c0 = (2.0 - b)*secderiv1;
0059   const double c1 = (1.0 + b)*secderiv2;
0060   return y1 + b*(y2 - y1) + (b*(b-1.0))*(c0+c1)*(dl*dl*os);
0061 }
0062
0063 double GetLinear(double x1, double x2, double y1, double y2, double x) {
0064   // Unchecked precondition: x1 < x < x2
0065   const double dl = x2 - x1;
0066   // note: all corner cases of the previous methods are covered and eventually
0067   //       gives b=0/1 that results in y=y0\y_{N-1} if e<=x[0]/e>=x[N-1] or
0068   //       y=y_i/y_{i+1} if e<x[i]/e>=x[i+1] due to small numerical errors
0069   const double  b = G4HepEmMax(0., G4HepEmMin(1., (x - x1)/dl));
0070   return y1 + b*(y2 - y1);
0071 }
0072
0073
0074 // use the improved, robust spline interpolation that I put in G4 10.6
0075 double GetSplineLog(int ndata, double* xdata, double* ydata, double* secderiv, double x, double logx, double logxmin, double invLDBin) {
0076   // make sure that $x \in  [x[0],x[ndata-1]]$
0077   const double xv = G4HepEmMax(xdata[0], G4HepEmMin(xdata[ndata-1], x));
0078   // compute the lowerindex of the x bin (idx \in [0,N-2] will be guaranted)
0079   const int   idx = (int)G4HepEmMax(0., G4HepEmMin((logx-logxmin)*invLDBin, ndata-2.));
0080   return GetSpline(xdata[idx], xdata[idx+1], ydata[idx], ydata[idx+1], secderiv[idx], secderiv[idx+1], xv);
0081 }
0082
0083 // same as above but both ydata and secderiv are stored in ydata array
0084 double GetSplineLog(int ndata, double* xdata, double* ydata, double x, double logx, double logxmin, double invLDBin) {
0085   // make sure that $x \in  [x[0],x[ndata-1]]$
0086   const double xv = G4HepEmMax(xdata[0], G4HepEmMin(xdata[ndata-1], x));
0087   // compute the lowerindex of the x bin (idx \in [0,N-2] will be guaranted)
0088   const int   idx = (int)G4HepEmMax(0., G4HepEmMin((logx-logxmin)*invLDBin, ndata-2.));
0089   const int  idx2 = 2*idx;
0090   return GetSpline(xdata[idx], xdata[idx+1], ydata[idx2], ydata[idx2+2], ydata[idx2+1], ydata[idx2+3], xv);
0091 }
0092
0093
0094 // same as above but all xdata, ydata and secderiv are stored in data array
0095 double GetSplineLog(int ndata, double* data, double x, double logx, double logxmin, double invLDBin) {
0096   // make sure that $x \in  [x[0],x[ndata-1]]$
0097   const double xv = G4HepEmMax(data[0], G4HepEmMin(data[3*(ndata-1)], x));
0098   // compute the lowerindex of the x bin (idx \in [0,N-2] will be guaranted)
0099   const int   idx = (int)G4HepEmMax(0., G4HepEmMin((logx-logxmin)*invLDBin, ndata-2.));
0100   const int  idx3 = 3*idx;
0101   return GetSpline(data[idx3], data[idx3+3], data[idx3+1], data[idx3+4], data[idx3+2], data[idx3+5], xv);
0102 }
0103
0104
0105 // this is used for getting inverse-range on host
0106 double GetSpline(double* xdata, double* ydata, double* secderiv, double x, int idx, int step) {
0107   return GetSpline(xdata[step*idx], xdata[step*(idx+1)], ydata[idx], ydata[idx+1], secderiv[idx], secderiv[idx+1], x);
0108 }
0109
0110 // same as above but both ydata and secderiv are stored in ydata array
0111 double GetSpline(double* xdata, double* ydata, double x, int idx) {
0112   const int  idx2 = 2*idx;
0113   return GetSpline(xdata[idx], xdata[idx+1], ydata[idx2], ydata[idx2+2], ydata[idx2+1], ydata[idx2+3], x);
0114 }
0115
0116 // same as above but both xdata, ydata and secderiv are stored in data array
0117 double GetSpline(double* data, double x, int idx) {
0118   const int  idx3 = 3*idx;
0119   return GetSpline(data[idx3], data[idx3+3], data[idx3+1], data[idx3+4], data[idx3+2], data[idx3+5], x);
0120 }
0121
0122
0123
0124
0125 double GetLinearLog(int ndata, double* data, double x, double logx, double logxmin, double invLDBin) {
0126   // make sure that $x \in  [x[0],x[ndata-1]]$
0127   //const double xv = G4HepEmMax(data[0], G4HepEmMin(data[2*(ndata-1)], x));
0128   // compute the lowerindex of the x bin (idx \in [0,N-2] will be guaranted)
0129   const int   idx = (int)G4HepEmMax(0., G4HepEmMin((logx-logxmin)*invLDBin, ndata-2.));
0130   const int  idx2_0 = 2*idx;
0131   const int  idx2_1 = idx2_0+2;
0132   return GetLinear(data[idx2_0], data[idx2_1], data[idx2_0+1], data[idx2_1+1], x);
0133 }
0134
0135 // interpolates one data out of the 2 specifed by the `iwhich=1,2` parameter
0136 double GetLinearLog2(int ndata, double* data, double x, double logx, double logxmin, double invLDBin, int iwhich) {
0137   // make sure that $x \in  [x[0],x[ndata-1]]$
0138   //const double xv = G4HepEmMax(data[0], G4HepEmMin(data[3*(ndata-1)], x));
0139   // compute the lowerindex of the x bin (idx \in [0,N-2] will be guaranted)
0140   const int   idx = (int)G4HepEmMax(0., G4HepEmMin((logx-logxmin)*invLDBin, ndata-2.));
0141   const int  idx3_0 = 3*idx;
0142   const int  idx3_1 = idx3_0+3;
0143   return GetLinear(data[idx3_0], data[idx3_1], data[idx3_0+iwhich], data[idx3_1+iwhich], x);
0144 }
0145 // interpolates all the 2 data at once
0146 void GetLinearLog2(int ndata, double* data, double x, double logx, double logxmin, double invLDBin, double res[2]) {
0147   // make sure that $x \in  [x[0],x[ndata-1]]$
0148   //const double xv = G4HepEmMax(data[0], G4HepEmMin(data[3*(ndata-1)], x));
0149   // compute the lowerindex of the x bin (idx \in [0,N-2] will be guaranted)
0150   const int   idx = (int)G4HepEmMax(0., G4HepEmMin((logx-logxmin)*invLDBin, ndata-2.));
0151   const int  idx3_0 = 3*idx;
0152   const int  idx3_1 = idx3_0+3;
0153   //return GetLinear(data[idx3_0], data[idx3_1], data[idx3_0+iwhich], data[idx3_1+iwhich], x);
0154   const double dl = data[idx3_1] - data[idx3_0];
0155   const double  b = G4HepEmMax(0., G4HepEmMin(1., (x - data[idx3_0])/dl));
0156   for (int i=1; i<3; ++i) {
0157     const double y1 = data[idx3_0+i];
0158     const double y2 = data[idx3_1+i];
0159     res[i-1] = G4HepEmMax(0.0, y1 + b*(y2 - y1));
0160   }
0161 }
0162 /*
0163 double GetLinearLogN(int ndata, double* data, double x, double logx, double logxmin, double invLDBin, int iwhich, int shift) {
0164   // make sure that $x \in  [x[0],x[ndata-1]]$
0165   const double xv = G4HepEmMax(data[0], G4HepEmMin(data[shift*(ndata-1)], x));
0166   // compute the lowerindex of the x bin (idx \in [0,N-2] will be guaranted)
0167   const int   idx = (int)G4HepEmMax(0., G4HepEmMin((logx-logxmin)*invLDBin, ndata-2.));
0168   const int    i0 = shift*idx;
0169   const int    i1 = shift*(idx+1);
0170   return GetLinear(data[i0], data[i1], data[i0+iwhich], data[i1+iwhich], x);
0171 }
0172 */
0173
0174
0175
0176
0177 // interpolates one out of the 4 data at once, specififed by `iwhich =1,2,3 or 4`
0178 double GetSplineLog4(int ndata, double* data, double x, double logx, double logxmin, double invLDBin, int iwhich) {
0179   // make sure that $x \in  [x[0],x[ndata-1]]$
0180   //const double xv = G4HepEmMax(data[0], G4HepEmMin(data[9*(ndata-1)], x));
0181   // compute the lowerindex of the x bin (idx \in [0,N-2] will be guaranted)
0182   const int   idx = (int)G4HepEmMax(0., G4HepEmMin((logx-logxmin)*invLDBin, ndata-2.));
0183   const int  idx9_0 = 9*idx;
0184   const int  idx9_1 = idx9_0 + 9;
0185   iwhich = (iwhich-1)*2+1;
0186   return GetSpline(data[idx9_0], data[idx9_1], data[idx9_0+iwhich], data[idx9_1+iwhich], data[idx9_0+iwhich+1], data[idx9_1+iwhich+1], x);
0187 }
0188 // interpolates all the 4 data at once
0189 void GetSplineLog4(int ndata, double* data, double x, double logx, double logxmin, double invLDBin, double res[4]) {
0190   // make sure that $x \in  [x[0],x[ndata-1]]$
0191   //const double xv = G4HepEmMax(data[0], G4HepEmMin(data[9*(ndata-1)], x));
0192   // compute the lowerindex of the x bin (idx \in [0,N-2] will be guaranted)
0193   const int   idx = (int)G4HepEmMax(0., G4HepEmMin((logx-logxmin)*invLDBin, ndata-2.));
0194   const int  idx9_0 = 9*idx;
0195   const int  idx9_1 = idx9_0 + 9;
0196   const double   dl = data[idx9_1] - data[idx9_0];
0197   const double    b = G4HepEmMax(0., G4HepEmMin(1., (x - data[idx9_0])/dl));
0198   for (int i=0; i<4; i++) {
0199     const double os = 0.166666666667; // 1./6.
0200     const int ii = 2*i+1;
0201     const double secderiv1 = data[idx9_0+ii+1];
0202     const double secderiv2 = data[idx9_1+ii+1];
0203     const double c0 = (2.0 - b)*secderiv1;
0204     const double c1 = (1.0 + b)*secderiv2;
0205     const double y1 = data[idx9_0+ii];
0206     const double y2 = data[idx9_1+ii];
0207     res[i] = G4HepEmMax(0.0, y1 + b*(y2 - y1) + (b*(b-1.0))*(c0+c1)*(dl*dl*os));
0208   }
0209 }
0210 /*
0211 double GetSplineLogN(int ndata, double* data, double x, double logx, double logxmin, double invLDBin, int iwhich, int shift) {
0212   // make sure that $x \in  [x[0],x[ndata-1]]$
0213   const double xv = G4HepEmMax(data[0], G4HepEmMin(data[shift*(ndata-1)], x));
0214   // compute the lowerindex of the x bin (idx \in [0,N-2] will be guaranted)
0215   const int   idx = (int)G4HepEmMax(0., G4HepEmMin((logx-logxmin)*invLDBin, ndata-2.));
0216   const int    i0 = shift*idx;
0217   const int    i1 = shift*(idx+1);
0218   iwhich = (iwhich-1)*2+1;
0219   return GetSpline(data[i0], data[i1], data[i0+iwhich], data[i1+iwhich], data[i0+iwhich+1], data[i1+iwhich+1], x);
0220 }
0221 */
0222
0223
0224 // this is used to get index for inverse range on host
0225 // NOTE: it is assumed that x[0] <= x and x < x[step*(num-1)]
0226 // step: the delta with which   the x values are located in xdata (i.e. =1 by default)
0227 int    FindLowerBinIndex(double* xdata, int num, double x, int step) {
0228   // Perform a binary search to find the interval val is in
0229   int ml = -1;
0230   int mu = num-1;
0231   while (std::abs(mu-ml)>1) {
0232     int mav = 0.5*(ml+mu);
0233     if (x<xdata[step*mav]) {  mu = mav; }
0234     else                   {  ml = mav; }
0235   }
0236   return mu-1;
0237 }