Index: /trunk/Mars/msimreflector/MFresnelLens.cc =================================================================== --- /trunk/Mars/msimreflector/MFresnelLens.cc (revision 19637) +++ /trunk/Mars/msimreflector/MFresnelLens.cc (revision 19638) @@ -27,6 +27,12 @@ // MFresnelLens // -// For some details on definitions please refer to -// https://application.wiley-vch.de/berlin/journals/op/07-04/OP0704_S52_S55.pdf +// For some details on definitions please refer to +// https://application.wiley-vch.de/berlin/journals/op/07-04/OP0704_S52_S55.pdf +// +// The HAWC's Eye lens is an Orafol SC943 +// https://www.orafol.com/en/europe/products/optic-solutions/productlines#pl1 +// +// A good description on ray-tracing can be found here +// https://graphics.stanford.edu/courses/cs148-10-summer/docs/2006--degreve--reflection_refraction.pdf // ////////////////////////////////////////////////////////////////////////////// @@ -41,4 +47,6 @@ #include "MReflection.h" +#include "MMath.h" + #include "MLog.h" #include "MLogManip.h" @@ -48,15 +56,273 @@ using namespace std; +// ========================================================================== + +enum exception_t +{ + kValidRay = 0, + + kStrayUpgoing, + kOutsideRadius, + kNoSurfaceFound, + kStrayDowngoing, + kAbsorbed, + + kFoundSurfaceUnavailable, + + kInvalidOrigin, + kTransitionError, + + kEnter = 1000, + kLeave = 2000, +}; + +enum surface_t +{ + kPhotonHasLeft = 0, + + kEntrySurface, + kSlopeSurface, + kDraftSurface, + kExitSurface, + + kMaterial = 5, + + kNoSurface = 9 +}; + + +class raytrace_exception : public runtime_error +{ +protected: + int fError; + int fOrigin; + int fSurface; + +public: + raytrace_exception(const int &_id, const int &_origin, const int &_surface, const string& what_arg) : + runtime_error(what_arg), fError(_id), fOrigin(_origin), fSurface(_surface) + { + } + + raytrace_exception(const int &_id, const int &_origin, const int &_surface, const char* what_arg) : + runtime_error(what_arg), fError(_id), fOrigin(_origin), fSurface(_surface) + { + } + + int id() const { return fError + fSurface*10 + fOrigin*100; } + int error() const { return fError; } + int origin() const { return fOrigin; } + int surface() const { return fSurface; } +}; + +class raytrace_error : public raytrace_exception +{ +public: + raytrace_error(const int &_id, const int &_origin, const int &_surface, const string& what_arg) : + raytrace_exception(_id, _origin, _surface, what_arg) { } + raytrace_error(const int &_id, const int &_origin, const int &_surface, const char* what_arg) : + raytrace_exception(_id, _origin, _surface, what_arg) { } +}; +class raytrace_info : public raytrace_exception +{ +public: + raytrace_info(const int &_id, const int &_origin, const int &_surface, const string& what_arg) : + raytrace_exception(_id, _origin, _surface, what_arg) { } + raytrace_info(const int &_id, const int &_origin, const int &_surface, const char* what_arg) : + raytrace_exception(_id, _origin, _surface, what_arg) { } +}; +class raytrace_user : public raytrace_exception +{ +public: + raytrace_user(const int &_id, const int &_origin, const int &_surface, const string& what_arg) : + raytrace_exception(_id, _origin, _surface, what_arg) { } + raytrace_user(const int &_id, const int &_origin, const int &_surface, const char* what_arg) : + raytrace_exception(_id, _origin, _surface, what_arg) { } +}; + +// ========================================================================== + + // -------------------------------------------------------------------------- // // Default constructor // -MFresnelLens::MFresnelLens(const char *name, const char *title) +MFresnelLens::MFresnelLens(const char *name, const char *title) : + fPSF(0), fSlopeAbsorption(false), fDraftAbsorption(false), + fBottomReflection(true), fDisableMultiEntry(false), fFresnelReflection(true), + fMinHits(0), fMaxHits(0) { fName = name ? name : "MFresnelLens"; fTitle = title ? title : "Parameter container storing a collection of several mirrors (reflector)"; - fMaxR = 55/2.; -} + // Default: Orafol SC943 + + DefineLens(); +} + +// ========================================================================== + +// -------------------------------------------------------------------------- +// +// Default ORAFOL SC943 +// +// Focal Length: F = 50.21 cm +// Diameter: D = 54.92 cm +// Groove width: w = 0.01 cm +// Lens thickness: h = 0.25 cm +// +// Default wavelength: 546 nm +// +void MFresnelLens::DefineLens(double F, double D, double w, double h, double lambda) +{ + fR = D/2; // [cm] Lens radius + fW = w; // [cm] Width of a single groove + fH = h; // [cm] Thickness of lens + fF = F; // [cm] focal length (see also MGeomCamFAMOUS!) + + fLambda = lambda; + + fN = MFresnelLens::RefractiveIndex(fLambda); // Lens + + // Velocity of light within the lens material [cm/ns] + // FIXME: Note that for the correct conversion in Transmission() + // also the speed in the surrounding medium has to be taken correctly + // into account (here it is assumed to be air with N=1 + fVc = fN/(TMath::C()*100/1e9); // cm/ns + + InitGeometry(fR, fW, fN, fF, fH); +} + +// -------------------------------------------------------------------------- +// +// Precalculate values such as the intersection points inside the grooves, +// the angle of the slope and draft surface and the corresponding tangents. +// +void MFresnelLens::InitGeometry(double maxr, double width, double N0, double F, double d) +{ + const uint32_t num = TMath::CeilNint(maxr/width); + + fGrooves.resize(num); + + for (int i=0; i. +// +// If no valid data was read, 0 is returned. -1 is returned if any tranmission +// value read from the file is >1. If the fresnel correction leads to a value >1, +// the value is set to 1. The number of valid data points is returned. +// +Int_t MFresnelLens::ReadTransmission(const TString &file, float thickness, bool correction) +{ + TGraph transmission(file); + + /* + double gx_min, gx_max, gy_min, gy_max; + absorption.ComputeRange(gx_min, gy_min, gx_max, gy_max); + if (lambdagx_max) + { + cout << "Invalid wavelength" << endl; + return; + }*/ + + if (transmission.GetN()==0) + return 0; + + for (int i=0; i1) + { + gLog << err << "Transmission larger than 1." << endl; + return -1; + } + + if (correction) + { + // Something like this is requried if correction + // for optical boundaries is necessary + const double n0 = MFresnelLens::RefractiveIndex(lambda); + + // FIXME: Make N_air a variable + const double r0 = (n0-1.0003)/(n0+1.0003); + const double r2 = r0*r0; + + trans *= (1+r2)*(1+r2); + + if (trans>1) + { + gLog << warn << "Transmission at " << lambda << "nm (" << trans << ") after Fresnel correction larger than 1." << endl; + trans = 1; + } + } + + // convert to absorption length (FIMXE: Sanity check) + transmission.GetY()[i] = -thickness/log(trans>0.999 ? 0.999 : trans); + } + + fAbsorptionLength = MSpline3(transmission); + + return fAbsorptionLength.GetNp(); +} + +// ========================================================================== // -------------------------------------------------------------------------- @@ -78,4 +344,58 @@ return sqrt(1+c0+c1+c2); +} + +// -------------------------------------------------------------------------- +// +// A Fresnel lens with parabolic surface calculated with the sagittag +// function (k=-1) and a correction for the thickness of the lens +// on the curvature. See also PhD thesis, Tim Niggemann ch. 7.1.1. +// +// Parameters are: +// The distance from the center r +// The focal length to be achieved F +// The refractive index of the outer medium (usually air) n0 +// The refractive index of the lens material (e.g. PMMA) n1 +// The thichness of the lens d +// +// r, F and d have to be in the same units. +// +// Return the slope angle alpha [rad]. The Slope angle is defined with +// respect to the plane of the lens. (0 at the center, decreasing +// with increasing radial distance) +// +double MFresnelLens::SlopeAngleParabolic(double r, double F, double n0, double n1, double d) +{ + // PhD, Tim Niggemann, ch 7.1.1 + const double rn = n1/n0; + const double c = (rn - 1) * (F + d/rn); + return -atan(r/c); +} + +// -------------------------------------------------------------------------- +// +// A Fresnel lens with an optimized parabolic surface calculated with +// the sagittag function (k=-1) and fitted coefficients according +// to Master thesis, Eichler. +// +// Note that for this setup other parameters must be fixed +// +// Parameters are: +// The distance from the center r +// +// r is in cm. +// +// Return the slope angle alpha [rad]. The Slope angle is defined with +// respect to the plane of the lens. (0 at the center, decreasing +// with increasing radial distance) +// +double MFresnelLens::SlopeAngleAspherical(double r) +{ + // Master, Eichler [r/cm] + return -atan( r/26.47 + +2*1.18e-4 * 1e1*r + +4*1.34e-9 * 1e3*r*r*r + +6*9.52e-15 * 1e5*r*r*r*r*r + -8*2.04e-19 * 1e7*r*r*r*r*r*r*r); } @@ -92,7 +412,4 @@ // Here, a thin lens is assumed // -// The HAWC's Eye lens is an Orafol SC943 -// https://www.orafol.com/en/europe/products/optic-solutions/productlines#pl1 -// // sin(omega) = r / sqrt(r^2+F^2) // tan(alpha) = sin(omega) / [ 1 - sqrt(n^2-sin(omega)^2) ] @@ -104,8 +421,18 @@ // distance) // -double MFresnelLens::SlopeAngle(double r, double F, double n) -{ +double MFresnelLens::SlopeAngleOptimized(double r, double F, double n) +{ + // Use F+d/2 double so = r / sqrt(r*r + F*F); return atan(so / (1-sqrt(n*n - so*so))); // alpha<0, Range [0deg; -50deg] +} + +// -------------------------------------------------------------------------- +// +// Currently calles SlopeAngleParabolic(r, F, 1, n, d) +// +double MFresnelLens::SlopeAngle(double r, double F, double n, double d) +{ + return SlopeAngleParabolic(r, F, 1.0003, n, d); } @@ -126,4 +453,6 @@ return (3 + r*0.473)*TMath::DegToRad(); // Range [0deg; 15deg] } + +// ========================================================================== // -------------------------------------------------------------------------- @@ -149,45 +478,57 @@ } -struct Groove -{ - double fTheta; - double fPeakZ; - - double fTanTheta; - double fTanTheta2[3]; - - Groove() { } - Groove(double theta, double z) : fTheta(theta), fPeakZ(z) - { - fTanTheta = tan(theta); - - fTanTheta2[0] = fTanTheta*fTanTheta; - fTanTheta2[1] = fTanTheta*fTanTheta * fPeakZ; - fTanTheta2[2] = fTanTheta*fTanTheta * fPeakZ*fPeakZ; - } -}; - -double CalcIntersection(const MQuaternion &p, const MQuaternion &u, const Groove &g, const double &fThickness) -{ -/* - Z: position of peak - - r_cone(z) = (Z-z)*tan(theta) - - (x) (p.x) (u.x) - (y) = (p.y) + dz * (u.y) - (z) (p.z) (u.z) - - r_line(z) = sqrt(x^2 + y^2) = sqrt( (p.x+dz*u.x)^2 + (p.y+dz*u.y)^2) - - Solved with wmmaxima: - - solve((H-z)^2*t^2 = (px+(z-pz)/uz*ux)^2+(py+(z-pz)/uz*uy)^2, z); - - [ - z=-(uz*sqrt((py^2+px^2)*t^2*uz^2+((2*H*py-2*py*pz)*t^2*uy+(2*H*px-2*px*pz)*t^2*ux)*uz+((pz^2-2*H*pz+H^2)*t^2-px^2)*uy^2+2*px*py*ux*uy+((pz^2-2*H*pz+H^2)*t^2-py^2)*ux^2)-H*t^2*uz^2+(-py*uy-px*ux)*uz+pz*uy^2+pz*ux^2)/(t^2*uz^2-uy^2-ux^2) - z= (uz*sqrt((py^2+px^2)*t^2*uz^2+((2*H*py-2*py*pz)*t^2*uy+(2*H*px-2*px*pz)*t^2*ux)*uz+((pz^2-2*H*pz+H^2)*t^2-px^2)*uy^2+2*px*py*ux*uy+((pz^2-2*H*pz+H^2)*t^2-py^2)*ux^2)+H*t^2*uz^2+( py*uy+px*ux)*uz-pz*uy^2-pz*ux^2)/(t^2*uz^2-uy^2-ux^2) - ] -*/ +// ========================================================================== + +// FIXME: The rays could be 'reflected' inside the material +// (even though its going out) or vice versa +static double RandomTheta(double psf) +{ + return psf>0 ? MMath::RndmPSF(psf)/2 : 0; +} + +// FIXME: The rays could be 'reflected' inside the material +// (even though its going out) or vice versa +static double RandomPhi(double r, double psf) +{ + return psf>0 ? MMath::RndmPSF(psf)/2 : 0; +} + + +// -------------------------------------------------------------------------- +// +// Calculate the intersection point beweteen a line defined by the position p +// and the direction u and a cone defined by the object cone. +// +// Z: position of peak of cone +// theta: opening angle of cone +// +// Distance r of cone surface at given z from z-axis +// r_cone(z) = (Z-z)*tan(theta) +// +// Equalition of line +// (x) (p.x) (u.x/u.z) +// (y) = (p.y) + dz * (u.y/u.z) +// (z) (p.z) ( 1 ) +// +// Normalization +// U.x := u.x/u.z +// U.y := u.y/u.z +// +// Distance of line at given z from z-axis +// r_line(z) = sqrt(x^2 + y^2) = sqrt( (p.x+dz*u.x)^2 + (p.y+dz*u.y)^2) with dz = z-p.z +// +// Equation to be solved +// r_cone(z) = r_line(z) +// +// Solved with wxmaxima: +// +// [0] solve((px+(z-pz)*Ux)^2+(py+(z-pz)*Uy)^2= ((Z-z)*t)^2, z); +// +// z= (sqrt(((Uy^2+Ux^2)*pz^2+(-2*Uy*py-2*Ux*px-2*Z*Uy^2-2*Z*Ux^2)*pz+py^2+2*Z*Uy*py+px^2+2*Z*Ux*px+Z^2*Uy^2+Z^2*Ux^2)*t^2-Ux^2*py^2+2*Ux*Uy*px*py-Uy^2*px^2)+Z*t^2+(-Uy^2-Ux^2)*pz+Uy*py+Ux*px)/(t^2-Uy^2-Ux^2), +// z=-(sqrt(((Uy^2+Ux^2)*pz^2+(-2*Uy*py-2*Ux*px-2*Z*Uy^2-2*Z*Ux^2)*pz+py^2+2*Z*Uy*py+px^2+2*Z*Ux*px+Z^2*Uy^2+Z^2*Ux^2)*t^2-Ux^2*py^2+2*Ux*Uy*px*py-Uy^2*px^2)-Z*t^2+( Uy^2+Ux^2)*pz-Uy*py-Ux*px)/(t^2-Uy^2-Ux^2) +// +double MFresnelLens::CalcIntersection(const MQuaternion &p, const MQuaternion &u, const Cone &cone) const +{ + const double &Z = cone.z; const double Ux = u.X()/u.Z(); @@ -198,8 +539,6 @@ const double pz = p.Z(); - const double H = g.fPeakZ; - - const double t = g.fTanTheta; - const double t2 = t*t; + //const double &t = cone.tan_theta; + const double &t2 = cone.tan_theta2; const double Ur2 = Ux*Ux + Uy*Uy; @@ -209,242 +548,1011 @@ const double cr2 = Ux*py - Uy*px; - const double a = t2 - Ur2; - const double b = Ur2*pz - Up2 - H*t2; - - const double h = H-pz; - const double h2 = h*h; + const double a = t2 - Ur2; + const double b = Ur2*pz - Up2 - Z*t2; + + const double h = Z-pz; + const double h2 = h*h; // [ -b +-sqrt(b^2 - 4 ac) ] / [ 2a ] const double radix = (Ur2*h2 + 2*Up2*h + pr2)*t2 - cr2*cr2; - if (radix<0) return 0; - double dz[2] = - { - (-b+sqrt(radix))/a, - (-b-sqrt(radix))/a + const double sqrt_radix = sqrt(radix); + + const double dz[2] = + { + (+sqrt_radix - b)/a, + (-sqrt_radix - b)/a }; - if (dz[0]<0 && dz[0]>fThickness) - return dz[0]; - - if (dz[1]<0 && dz[1]>fThickness) - return dz[1]; - - return 0; + // Return the closest solution inside the allowed range + // which is in the direction of movement + + const double &H = cone.h; + + const bool is_inside0 = dz[0]>=H && dz[0]<0; + const bool is_inside1 = dz[1]>=H && dz[1]<0; + + // FIXME: Simplify! + if (!is_inside0 && !is_inside1) + return 0; + + // Only dz[0] is in the right z-range + if (is_inside0 && !is_inside1) + { + // Check if dz[0] is in the right direction + if ((u.Z()>=0 && dz[0]>=p.Z()) || + (u.Z()< 0 && dz[0]< p.Z())) + return dz[0]; + + return 0; + } + + // Only dz[1] is in the right z-range + if (!is_inside0 && is_inside1) + { + // Check if dz[1] is in the right direction + if ((u.Z()>=0 && dz[1]>=p.Z()) || + (u.Z()< 0 && dz[1]< p.Z())) + return dz[1]; + + return 0; + } + + /* + if (is_inside0^is_inside1) + { + if (u.Z()>=0) + return dz[0]>p.Z() ? dz[0] : dz[1]; + else + return dz[0]=0) + { + // Both solution could be correct + if (dz[0]>=p.Z() && dz[1]>=p.Z()) + return std::min(dz[0], dz[1]); + + // only one solution can be correct + return dz[0]>=p.Z() ? dz[0] : dz[1]; + } + else + { + // Both solution could be correct + if (dz[0] 0; // r is (at least locally) increasing + + // --------------------------- + const double Ux = u.X()/u.Z(); + const double Uy = u.Y()/u.Z(); + + const double px = p.X(); + const double py = p.Y(); + const double pz = p.Z(); + + const double Ur2 = Ux*Ux + Uy*Uy; + const double cr2 = Ux*py - Uy*px; + const double pr2 = px*px + py*py; + const double Up2 = Ux*px + Uy*py; + + //for (int i=1; i=z2 && dz[0]<=z1) + return i; + + if (dz[1]>=z2 && dz[1]<=z1) + return i; + } + } + + return -1; +} + +// -------------------------------------------------------------------------- +// +// If no transmission was given returns true. Otherwaise calculates the +// absorption length for a flight time dt in the material and a photon +// with wavelength lambda. The flight time is converted to a geometrical +// using the speed of light in the medium. +// +// Returns true if the poton passed, false if it was absorbed. +// +bool MFresnelLens::Transmission(double dt, double lambda) const +{ + if (fAbsorptionLength.GetNp()==0) + return true; + + // FIXME: Speed up! + const double alpha = fAbsorptionLength.Eval(lambda); + + // We only have the travel time, thus we have to convert back to distance + // Note that the transmission coefficients are w.r.t. to geometrical + // distance not light-travel distance. Thus the distance has to be corrected + // for the corresponding refractive index of the material. + const double cm = dt/fVc; + + const double trans = exp(-cm/alpha); + return gRandom->Uniform()=fR) + { + //cout << "Left the lens radius (enter)" << endl; + throw -2; + } + //if (dir.Z()>0) + //{ + // cout << "Upgoing, outside of the material" << endl; + // PropagateZ(pos, dir, dir.Z()>0 ? 3 : -3); + // return -1; + //} + + + // Calculate the ordinal number of the groove correpsonding to rx + const int ix = TMath::FloorNint(rx/fW); + + // Photons was just injected (test both surfaces) or came from the other surface + if (surface==0 || surface==2) + { + // Get the possible intersection point with the slope angle + const double z1 = CalcIntersection(pos, dir, fGrooves[ix].slope); + + // We hit the slope angle + if (z1!=0) + { + // Move photon to new hit position + pos.PropagateZ(dir, z1); + + if (fSlopeAbsorption) + throw -100; + + // Get the normal vector of the surface which was hit + const VectorNorm norm(fGrooves[ix].slope.theta_norm+RandomTheta(fPSF), + pos.XYvector().Phi()+RandomPhi(pos.R(), fPSF)); + + // Get the optical transition of the direction vector + const int ret = MOptics::ApplyTransition(dir, norm, 1, n0); + + // Transition was Reflection - try again + if (ret==1 || ret==2) + return EnterGroove(1, n0, lambda, pos, dir)+1; + + // Transition was Refraction - enter + if (ret>=3) + return LeavePeak(1, n0, lambda, pos, dir, pos.T())+1; + + // Error occured (see ApplyTransition for details) + //cout << "ERR[TIR1]" << endl; + throw -3; + } + } + + // Photons was just injected (test both surfaces) or came from the other surface + if (surface==0 || surface==1) + { + const double z2 = CalcIntersection(pos, dir, fGrooves[ix].draft); + + // We hit the draft angle + if (z2!=0) + { + // Move photon to new hit position + pos.PropagateZ(dir, z2); + + if (fDraftAbsorption) + throw -101; + + // Get the normal vector of the surface which was hit + const VectorNorm norm(fGrooves[ix].draft.theta_norm+RandomTheta(fPSF), + pos.XYvector().Phi()+RandomPhi(pos.R(), fPSF)); + + // Get the optical transition of the direction vector + const int ret = MOptics::ApplyTransition(dir, norm, 1, n0); + + // Transition was Reflection - try again + if (ret==1 || ret==2) + return EnterGroove(2, n0, lambda, pos, dir)+1; + + // Transition was Refraction - enter + if (ret>=3) + return -LeavePeak(2, n0, lambda, pos, dir, pos.T())+1; + + // Error occured (see ApplyTransition for details) + //cout << "ERR[TIR2]" << endl; + throw -4; + } + } + + if (dir.Z()>0) + { + //cout << "Upgoing, outside of the material" << endl; + //pos.PropagateZ(dir, dir.Z()>0 ? 3 : -3); + throw -5; + } + + // The ray has left the peak at the bottom(?) + //cout << "ERR[N/A]" << endl; + throw -6; +} +*/ + + +// surface=0 : incoming ray +// surface=1 : slope +// surface=2 : draft +// surface=3 : bottom +int MFresnelLens::EnterGroove(int surface, double n0, MQuaternion &pos, MQuaternion &dir) const +{ + const double rx = pos.R(); + + if (surface==kExitSurface) + throw raytrace_error(kEnter+kInvalidOrigin, surface, -1, + "EnterGroove - Bottom as origin invalid"); + + if (rx>=fR) // This is an error as the direction vector is now invalid + throw raytrace_error(kEnter+kOutsideRadius, surface, -1, + "EnterGroove - Surface hit outside allowed radius"); + + /* + if (dir.Z()>0) + return -1; + }*/ + + + // FIXME: There is a very tiny chance that a ray hits the same surface twice for + // very horizontal rays. Checking this needs to make sure that the same + // solution is not just found again. + + // Calculate the ordinal number of the groove correpsonding to rx + const int ix = TMath::FloorNint(rx/fW); + + // Photons was just injected (test both surfaces) or came from the other surface + if (surface==kEntrySurface || surface==kDraftSurface) + { + // Get the possible intersection point with the slope angle + const double z1 = CalcIntersection(pos, dir, fGrooves[ix].slope); + + // We hit the slope angle + if (z1!=0) + { + // Move photon to new hit position + pos.PropagateZ(dir, z1); + if (fSlopeAbsorption) + throw raytrace_user(kEnter+kAbsorbed, surface, kSlopeSurface, + "EnterGroove - Photon absorbed by slope surface"); + + // Get the normal vector of the surface which was hit + const VectorNorm norm(fGrooves[ix].slope.theta_norm+RandomTheta(fPSF), + pos.XYvector().Phi()+RandomPhi(pos.R(), fPSF)); + + // Get the optical transition of the direction vector + const int ret = MOptics::ApplyTransition(dir, norm, 1, n0, fFresnelReflection); + + // Transition was Reflection - try again + if (ret==1 || ret==2) + return kSlopeSurface;//EnterGroove(1, n0, lambda, pos, dir)+1; + + // Transition was Refraction - enter + if (ret>=3) + return -kSlopeSurface;//LeavePeak(1, n0, lambda, pos, dir, pos.T())+1; + + // Error occured (see ApplyTransition for details) + throw raytrace_error(kEnter+kTransitionError, surface, kSlopeSurface, + "EnterGroove - MOptics::ApplyTransition failed for slope surface"); + } + } + + // Photons was just injected (test both surfaces) or came from the other surface + if (surface==kEntrySurface || surface==kSlopeSurface) + { + const double z2 = CalcIntersection(pos, dir, fGrooves[ix].draft); + + // We hit the draft angle + if (z2!=0) + { + // Move photon to new hit position + pos.PropagateZ(dir, z2); + if (fDraftAbsorption) + throw raytrace_user(kEnter+kAbsorbed, surface, kDraftSurface, + "EnterGroove - Photon absorbed by draft surface"); + + // Get the normal vector of the surface which was hit + const VectorNorm norm(fGrooves[ix].draft.theta_norm+RandomTheta(fPSF), + pos.XYvector().Phi()+RandomPhi(pos.R(), fPSF)); + + // Get the optical transition of the direction vector + const int ret = MOptics::ApplyTransition(dir, norm, 1, n0, fFresnelReflection); + + // Transition was Reflection - try again + if (ret==1 || ret==2) + return kDraftSurface;//EnterGroove(2, n0, lambda, pos, dir)+1; + + // Transition was Refraction - enter + if (ret>=3) + return -kDraftSurface;//LeavePeak(2, n0, lambda, pos, dir, pos.T())+1; + + // Error occured (see ApplyTransition for details) + throw raytrace_error(kEnter+kTransitionError, surface, kDraftSurface, + "EnterGroove - MOptics::ApplyTransition failed for draft surface"); + } + } + + if (dir.Z()>0) + { + // We have missed both surfaces and we are upgoing... + // ... ray can be discarded + throw raytrace_info(kEnter+kStrayUpgoing, surface, kNoSurface, + "EnterGroove - Particle is upgoing and has hit no surface"); + } + + // The ray has left the peak at the bottom(?) + throw raytrace_error(kEnter+kStrayDowngoing, surface, kNoSurface, + "EnterGroove - Particle is downgoing and has hit no surface"); +} + +/* +// Leave the peak from inside the material, either thought the draft surface or the +// slope surface or the bottom connecting the valley of both +int MFresnelLens::LeavePeak(int surface, double n0, double lambda, MQuaternion &pos, MQuaternion &dir, double T0) const +{ + const double rx = pos.R(); + + if (rx>=fR) + { + //cout << "Left the lens radius (leave)" << endl; + throw -10; + } + + if (dir.Z()>0 && surface!=3) // && surface!=4) + { + //cout << "Upgoing, inside of the material" << endl; + //pos.PropagateZ(dir, dir.Z()>0 ? 3 : -3); + throw -11; + } + + if (surface!=1 && surface!=2 && surface!=3) // && surface!=4) + { + //cout << "Surface of origin invalid" << endl; + throw -12; + } + + + // Calculate the ordinal number of the groove correpsonding to rx + const int ix = TMath::FloorNint(rx/fW); + + // FIXME: The Z-coordinate (cone.h) is actually a line through two points!!! + + Cone slope = fGrooves[ix].slope; + Cone draft = fGrooves[ix].draft; + + const bool is_draft = rx>fGrooves[ix].r; + if (is_draft) + { + // We are in the volume under the draft angle... taking the slope from ix+1 + if (ix0) // FIXME: Check whether this is correct + draft = fGrooves[ix-1].draft; + } + + if (is_draft+1!=surface && (surface==1 || surface==2)) + cout << "SURFACE: " << is_draft+1 << " " << surface << endl; + + if (surface==3) + { + //cout << "Upgoing, coming from the bottom of the lens" << endl; + // Find out which triangle (peak) the photon is going to enter + // then proceed... + throw -13; + } + + + // We are inside the material and downgoing, so if we come from a slope surface, + // we can only hit a draft surface after and vice versa + if (is_draft || surface==3) + { + const double z1 = CalcIntersection(pos, dir, slope); + + // We hit the slope angle and are currently in the volume under the draft surface + if (z1!=0) + { + // Move photon to new hit position + pos.PropagateZ(dir, z1); + + if (fSlopeAbsorption) + throw -200; + + // Get the normal vector of the surface which was hit + const VectorNorm norm(slope.theta_norm+RandomTheta(fPSF), + pos.XYvector().Phi()+RandomPhi(pos.R(), fPSF)); + + // Get the optical transition of the direction vector + const int ret = MOptics::ApplyTransition(dir, norm, n0, 1); + + // Transition was Reflection - try again + if (ret==1 || ret==2) + return LeavePeak(1, n0, lambda, pos, dir, T0)+1; + + // Transition was Refraction - leave + if (ret>=3) + { + if (!Transmission(pos.T()-T0, lambda)) + throw -14; + + return EnterGroove(1, n0, lambda, pos, dir)+1; + } + + // Error occured (see ApplyTransition for details) + //cout << "ERR[TIR3]" << endl; + throw -15; + } + } + + if (!is_draft || surface==3) + { + const double z2 = CalcIntersection(pos, dir, draft); + + // We hit the draft angle from the inside and are currently in the volume under the slope angle + if (z2!=0) + { + // Move photon to new hit position + pos.PropagateZ(dir, z2); + + if (fDraftAbsorption) + throw -201; + + // Get the normal vector of the surface which was hit + const VectorNorm norm(draft.theta_norm+RandomTheta(fPSF), + pos.XYvector().Phi()+RandomPhi(pos.R(), fPSF)); + + // Get the optical transition of the direction vector + const int ret = MOptics::ApplyTransition(dir, norm, n0, 1); + + // Transition was Reflection - try again + if (ret==1 || ret==2) + return LeavePeak(2, n0, lambda, pos, dir, T0)+1; + + // Transition was Refraction - leave + if (ret>=3) + { + if (!Transmission(pos.T()-T0, lambda)) + throw -16; + + return EnterGroove(2, n0, lambda, pos, dir)+1; + } + + // Error occured (see ApplyTransition for details) + //cout << "ERR[TIR4]" << endl; + throw -17; + } + } + + if (surface==3)// || surface==4) + { + //cout << ix << " Lost bottom reflected ray " << surface << endl; + throw -18; + } + + // The ray has left the peak at the bottom + + // FIXME: There is a tiny chance to escape to the side + // As there is a slope in the bottom surface of the peak + + // Move photon to new hit position + pos.PropagateZ(dir, -fH); + + if (pos.R()>fR) + { + //cout << "Left the lens radius (bottom)" << endl; + throw -19; + } + + // Get the normal vector of the surface which was hit + const VectorNorm norm(RandomTheta(fPSF), gRandom->Uniform(0, TMath::TwoPi())); + + // Get the optical transition of the direction vector + const int ret = MOptics::ApplyTransition(dir, norm, n0, 1); + + // Transition was Reflection + // (Photon scattered back from the bottom of the lens) + if (ret==1 || ret==2) + return LeavePeak(3, n0, lambda, pos, dir, T0)+1; + + // Transition was Refraction + // (Photon left at the bottom of the lens) + if (ret>=3) + { + if (!Transmission(pos.T()-T0, lambda)) + throw -20; + return 0; - - const double sqrtBac = sqrt(B2 - a*c); - - if (a==0) - return 0; - - const int sign = g.fPeakZ>=0 ? 1 : -1; - - const double dz[2] = - { - (sqrtBac+B) / a * sign, - (sqrtBac-B) / a * sign - }; - - // Only one dz[i] can be within the lens (due to geometrical reasons - - if (dz[0]<0 && dz[0]>fThickness) - return dz[0]; - - if (dz[1]<0 && dz[1]>fThickness) - return dz[1]; - - return 0; -} + } + + // Error occured (see ApplyTransition for details) + //cout << "ERR[TIR5]" << endl; + throw -21; +}*/ + +// Leave the peak from inside the material, either thought the draft surface or the +// slope surface or the bottom connecting the valley of both +int MFresnelLens::LeavePeak(int surface, double n0, MQuaternion &pos, MQuaternion &dir, double T0) const +{ + const double rx = pos.R(); + + if (rx>=fR) // This is an error as the direction vector is now invalid + throw raytrace_error(kLeave+kOutsideRadius, surface, kNoSurface, + "LeavePeak - Surface hit outside allowed radius"); + + // FIXME: Can we track them further? + if (fDisableMultiEntry && dir.Z()>0 && surface!=3/* && surface!=4*/) + throw raytrace_info(kLeave+kStrayUpgoing, surface, kNoSurface, + "LeavePeak - Particle is upgoing inside the material and does not come from the bottom"); + + if (surface!=kSlopeSurface && surface!=kDraftSurface && surface!=kExitSurface/* && surface!=4*/) + throw raytrace_error(kLeave+kInvalidOrigin, surface, kNoSurface, + "LeavePeak - Invalid surface of origin"); + + + // Calculate the ordinal number of the groove correpsonding to rx + const int ix = TMath::FloorNint(rx/fW); + + // FIXME: The Z-coordinate (cone.h) is actually a line through two points!!! + + Cone slope = fGrooves[ix].slope; + Cone draft = fGrooves[ix].draft; + + //if (is_draft+1!=surface && (surface==1 || surface==2)) + // cout << "SURFACE: " << is_draft+1 << " " << surface << endl; + + const bool is_draft = rx>fGrooves[ix].r; + if (is_draft) + { + // We are in the volume under the draft angle... taking the slope from ix+1 + if (ix0) // FIXME: Check whether this is correct + draft = fGrooves[ix-1].draft; + } + + if (surface==kExitSurface) + { + if (!fBottomReflection) + throw raytrace_user(kLeave+kAbsorbed, surface, kExitSurface, + "LeavePeak - Particle absorbed on the bottom"); + + const int in = FindPeak(ix, pos, dir); + + // This might happen if the ray is very flat and leaving + // the lens before hitting the border boundary of the grooves + if (in<0) + throw raytrace_error(kLeave+kNoSurfaceFound, kExitSurface, kNoSurface, + "LeavePeak - No hit surface found for particle reflected at the bottom"); + + slope = fGrooves[in].slope; + draft = fGrooves[in==0 ? 0 : in-1].draft; + } + + // FIXME: There is a chance that we can hit the same surface twice (for very horizontal rays + // but this requires a proper selection of the hit point + + // We are inside the material and downgoing, so if we come from a slope surface, + // we can only hit a draft surface after and vice versa + if (is_draft || surface==kExitSurface) + { + const double z1 = CalcIntersection(pos, dir, slope); + + // We hit the slope angle and are currently in the volume under the draft surface + if (z1!=0) + { + // Move photon to new hit position + pos.PropagateZ(dir, z1); + + if (fSlopeAbsorption) + throw raytrace_user(kLeave+kAbsorbed, surface, kSlopeSurface, + "LeavePeak - Photon absorbed by slope surface"); + + // Get the normal vector of the surface which was hit + const VectorNorm norm(slope.theta_norm+RandomTheta(fPSF), + pos.XYvector().Phi()+RandomPhi(pos.R(), fPSF)); + + // Get the optical transition of the direction vector + const int ret = MOptics::ApplyTransition(dir, norm, n0, 1, fFresnelReflection); + + // Transition was Reflection - try again + if (ret==1 || ret==2) + return -kSlopeSurface;//LeavePeak(1, n0, lambda, pos, dir, T0)+1; + + // Transition was Refraction - leave + if (ret>=3) // Transmission + return kSlopeSurface;//EnterGroove(1, n0, lambda, pos, dir)+1; + + // Error occured (see ApplyTransition for details) + throw raytrace_error(kLeave+kTransitionError, surface, kSlopeSurface, + "LeavePeak - MOptics::ApplyTransition failed for slope surface"); + } + } + + if (!is_draft || surface==kExitSurface) + { + const double z2 = CalcIntersection(pos, dir, draft); + + // We hit the draft angle from the inside and are currently in the volume under the slope angle + if (z2!=0) + { + // Move photon to new hit position + pos.PropagateZ(dir, z2); + + if (fDraftAbsorption) + throw raytrace_user(kLeave+kAbsorbed, surface, kDraftSurface, + "LeavePeak - Photon absorbed by draft surface"); + + // Get the normal vector of the surface which was hit + const VectorNorm norm(draft.theta_norm+RandomTheta(fPSF), + pos.XYvector().Phi()+RandomPhi(pos.R(), fPSF)); + + // Get the optical transition of the direction vector + const int ret = MOptics::ApplyTransition(dir, norm, n0, 1, fFresnelReflection); + + // Transition was Reflection - try again + if (ret==1 || ret==2) + return -kDraftSurface;//LeavePeak(2, n0, lambda, pos, dir, T0)+1; + + // Transition was Refraction - leave + if (ret>=3) // Transmission + return kDraftSurface;//EnterGroove(2, n0, lambda, pos, dir)+1; + + // Error occured (see ApplyTransition for details) + //cout << "ERR[TIR4]" << endl; + throw raytrace_error(kLeave+kTransitionError, surface, kDraftSurface, + "LeavePeak - MOptics::ApplyTransition failed for draft surface"); + } + } + + if (surface==kExitSurface/* || surface==4*/) + throw raytrace_error(kLeave+kFoundSurfaceUnavailable, kExitSurface, is_draft?kSlopeSurface:kDraftSurface, + "LeavePeak - Ray reflected on the bottom did not hit the found surface"); + + // The ray has left the peak at the bottom + + // FIXME: There is a tiny chance to escape to the side + // As there is a slope in the bottom surface of the peak + + // FIXME: Theoretically, a ray can hit the same surface twice + + // Move photon to new hit position + pos.PropagateZ(dir, -fH); + + if (pos.R()>fR) + throw raytrace_info(kLeave+kOutsideRadius, surface, kExitSurface, + "LeavePeak - Hit point at the bottom surface is beyond allowed radius"); + + // Get the normal vector of the surface which was hit + const VectorNorm norm(RandomTheta(fPSF), gRandom->Uniform(0, TMath::TwoPi())); + + // Get the optical transition of the direction vector + const int ret = MOptics::ApplyTransition(dir, norm, n0, 1, fFresnelReflection); + + // Transition was Reflection + // (Photon scattered back from the bottom of the lens) + if (ret==1 || ret==2) + return -kExitSurface;//LeavePeak(3, n0, lambda, pos, dir, T0)+1; + + // Transition was Refraction + // (Photon left at the bottom of the lens) + if (ret>=3) // Transmission + return kPhotonHasLeft; + + // Error occured (see ApplyTransition for details) + throw raytrace_error(kLeave+kTransitionError, surface, kExitSurface, "LeavePeak - MOptics::ApplyTransition failed for bottom surface"); +} + + +// Differences: +// Returns a 'reflected' vector at z=0 +// Does not propagate to z=0 at the beginning +Int_t MFresnelLens::ExecuteOptics(MQuaternion &p, MQuaternion &u, const Short_t &wavelength) const +{ + // Corsika Coordinates are in cm! + + const double lambda = wavelength==0 ? fLambda : wavelength; + if (fAbsorptionLength.GetNp()!=0 && + (lambdafAbsorptionLength.GetXmax())) + { + gLog << err << "Wavelength " << lambda << "nm out of absorption range [" << fAbsorptionLength.GetXmin() << "nm;" << fAbsorptionLength.GetXmax() << "nm]" << endl; + return -1; + } + + const double n0 = MFresnelLens::RefractiveIndex(lambda); + + try + { + int last_surface = EnterGroove(kEntrySurface, n0, p, u); + + // last_surface that was hit (photon originates from) + // 0 entrance (Z=0) or exit (Z=-fH) surface + // 1 slope + // 2 draft + // 3 bottom + // positive: photon is outside of material --> Try to enter + // nagative: photon is inside of material --> Try to leave + + double T0 = 0; + + // The general assumption is: no surface can be hit twice in a row + + int cnt = 0; + while (last_surface!=0) + { + cnt ++; + + // photon is outside of material --> try to enter + if (last_surface>0) + { + last_surface = EnterGroove( last_surface, n0, p, u); + + // successfully entered --> remember time of entrance to calculate transimission + if (last_surface<0) + T0 = p.T(); + + continue; + } + + // photon is inside of material --> try to leave + if (last_surface<0) + { + last_surface = LeavePeak(-last_surface, n0, p, u, T0); + + // successfully left --> apply transmission + if (last_surface>=0) + { + if (!Transmission(p.T()-T0, lambda)) + throw raytrace_error(kAbsorbed, last_surface, kMaterial, + "TraceRay - Ray absorbed in material"); + } + + continue; + } + } + + // To make this consistent with a mirror system, + // we now change our coordinate system + // Rays from the lens to the camera are up-going (positive sign) + u.fVectorPart.SetZ(-u.Z()); + + // In the datasheet, it looks as if F is calculated + // towards the center of the lens + // (Propagating to F means not propagating a distance of F-H/2) + p.fVectorPart.SetZ(0); + + return cnt>=fMinHits && (fMaxHits==0 || cnt<=fMaxHits) ? cnt : -1;; + } + catch (const raytrace_exception &e) + { + return -e.id(); + } + +/* + try + { + const int cnt = EnterGroove(0, n0, lambda, p, u); + + // To make this consistent with a mirror system, + // we now change our coordinate system + // Rays from the lens to the camera are up-going (positive sign) + u.fVectorPart.SetZ(-u.Z()); + + // In the datasheet, it looks as if F is calculated + // towards the center of the lens + // (Propagating to F means not propagating a distance of F-H/2) + p.fVectorPart.SetZ(0); + + return cnt>=fMinHits && (fMaxHits==0 || cnt<=fMaxHits) ? cnt : -1; + + } + catch (const int &rc) + { + return rc; + } */ - -Int_t MFresnelLens::ExecuteOptics(MQuaternion &p, MQuaternion &u, const Short_t &wavelength) const +} + +// Differences: +// Does propagate to z=0 at the beginning +Int_t MFresnelLens::TraceRay(vector &vec, MQuaternion &p, MQuaternion &u, const Short_t &wavelength, bool verbose) const { // Corsika Coordinates are in cm! - const double D = 54.92; // [cm] Diameter - const double F = 50.21; // [cm] focal length (see also MGeomCamFAMOUS!) - - const double R = D/2; // [cm] radius of lens - const double w = 0.01; // [cm] Width of a single groove - const double H = 0.25; // [cm] Thickness (from pdf) - - // Focal length is given at this wavelength - const double n0 = RefractiveIndex(546); - const double n = RefractiveIndex(wavelength==0 ? 546 : wavelength); - - const double r = p.R(); - - // Ray has missed the lens - if (r>R) - return -1; - - // Calculate the ordinal number of the groove which is hit - const unsigned int nth = TMath::FloorNint(r/w); - - // Calculate its lower and upper radius and its center - const double r0 = nth*w; // Lower radius of nth groove - const double rc = nth*w + w/2; // Upper radius of nth groove - const double r1 = nth*w + w; // Center of nth groove - - // Calculate the slope and draft angles - const double alpha = -SlopeAngle(rc, F, n0); - const double psi = DraftAngle(r1); - - // Define the corresponding surfaces - const Groove slope(TMath::Pi()/2-alpha, r0*tan(alpha)); - const Groove draft(-psi, -r1/tan(psi)); - - // Calculate the insersection of the ray between the two surfaces - // with z between 0 and -H - const double dz_slope = CalcIntersection(p, u, slope, -H); - const double dz_draft = CalcIntersection(p, u, draft, -H); - - // valid means that the photon has hit the fresnel surface, - // otherwise the draft surface - bool valid = dz_slope>dz_draft || dz_draft>=0; - - // Theta angle of the normal vector of the surface that was hit - TVector3 norm; - - // Propagate to the hit point. To get the correct normal vector of - // the surface, the hit point must be knwone - p.PropagateZ(u, valid ? dz_slope : dz_draft); - norm.SetMagThetaPhi(1, valid ? alpha : psi-TMath::Pi()/2, p.fVectorPart.Phi()); - - // Estimate reflectivity for this particular hit (n1 check before) - // Probability that the ray gets reflected - if (gRandom->Uniform()0) - return valid ? -3 : -4; - - // Propgate back to z=0 (FIXME: CalcIntersection requires z=0) - p.PropagateZ(u, 0); - - // Calc new intersection the other of the two surfaces - valid = !valid; - - const double dz = CalcIntersection(p, u, valid ? slope : draft, -H); - - // Propagate to the hit point at z=dz (dz<0) and get the new normal vector - p.PropagateZ(u, dz); - norm.SetMagThetaPhi(1, valid ? alpha : psi-TMath::Pi()/2, p.fVectorPart.Phi()); - } - - const double check_r = p.R(); - - // Some sanity checks (we loose a few permille due to numerical uncertainties) - // If the ray has not hit within the right radii.. reject - if (check_rr1) - return valid ? -7 : -8; - - // Find dw value of common z-value - // - // gz*tan(psi) = w - gw - // gz = dw*tan(alpha) - // - const double Gw = w/(tan(alpha)*tan(psi)+1); - const double Gz = -Gw*tan(alpha); - - if (valid && check_r>r0+Gw) - return -9; - if (!valid && check_rn2 => check after) - // Probability that the ray gets reflected - // (total internal reflection -> we won't track that further) - if (gRandom->Uniform()fAbsorptionLength.GetXmax())) + { + gLog << err << "Wavelength " << lambda << "nm out of absorption range [" << fAbsorptionLength.GetXmin() << "nm;" << fAbsorptionLength.GetXmax() << "nm]" << endl; + return -1; + } + + const double n0 = MFresnelLens::RefractiveIndex(lambda); + + // Photon must be at the lens surface + p.PropagateZ(u, 0); + vec.push_back(p); + + try + { + int last_surface = EnterGroove(kEntrySurface, n0, p, u); + //cout << "enter1 = " << last_surface << endl; + + // last_surface that was hit (photon originates from) + // 0 entrance (Z=0) or exit (Z=-fH) surface + // 1 slope + // 2 draft + // 3 bottom + // positive: photon is outside of material --> Try to enter + // nagative: photon is inside of material --> Try to leave + + double T0 = 0; + + // The general assumption is: no surface can be hit twice in a row + + int cnt = 0; + while (last_surface!=0) + { + cnt ++; + vec.push_back(p); + + // photon is outside of material --> try to enter + if (last_surface>0) + { + last_surface = EnterGroove( last_surface, n0, p, u); + //cout << "enter = " << last_surface << endl; + + // successfully entered --> remember time of entrance to calculate transimission + if (last_surface<0) + T0 = p.T(); + + continue; + } + + // photon is inside of material --> try to leave + if (last_surface<0) + { + last_surface = LeavePeak(-last_surface, n0, p, u, T0); + //cout << "leave = " << last_surface << endl; + + // successfully left --> apply transmission + if (last_surface>=0) + { + if (!Transmission(p.T()-T0, lambda)) + throw raytrace_error(kAbsorbed, last_surface, kMaterial, + "TraceRay - Ray absorbed in material"); + } + + continue; + } + } + + vec.push_back(p); + return cnt; + } + catch (const raytrace_exception &e) + { + if (verbose) + gLog << all << e.id() << ": " << e.what() << endl; + + // Hit point at bottom surface beyond allowed range + // FIXME: Only if surface is kExitSurface + if (e.id()==2342) + vec.push_back(p); + + return -e.id(); + } +} Index: /trunk/Mars/msimreflector/MFresnelLens.h =================================================================== --- /trunk/Mars/msimreflector/MFresnelLens.h (revision 19637) +++ /trunk/Mars/msimreflector/MFresnelLens.h (revision 19638) @@ -1,18 +1,103 @@ #ifndef MARS_MFresnelLens #define MARS_MFresnelLens + +#include #ifndef MARS_MOptics #include "MOptics.h" #endif - -class TVector3; -class MQuaternion; +#ifndef MARS_MSpline3 +#include "MSpline3.h" +#endif +#ifndef MARS_MQuaternion +#include "MQuaternion.h" +#endif class MFresnelLens : public MOptics { +public: + // --------------- Helper class --------------------- + // This is a TVector3 which is nromalized and takes theta and phi as arguent + // same as calling + // TVector3 v; + // v.SetMagThetaPhi(1, theta, phi); + class VectorNorm : public TVector3 + { + public: + VectorNorm(double theta, double phi) : + TVector3( + sin(theta) * cos(phi), + sin(theta) * sin(phi), + cos(theta) + ) + { + } + }; + + // ----------- Properties of the grooves ------------- + + struct Cone + { + double z; // z-coordinate of the peak of the cone (x=y=0) + double theta; // Opening angle of the cone [rad] + double tan_theta; // tan(theta) + double tan_theta2; // tan(theta)*tan(theta) + double h; // height of the slope (intersection point in the valley), h<0 + double theta_norm; // theta angle of the normal vector corresponding to the cone surface + }; + + struct Groove + { + double r; // Radius of the intersection point between slope and draft angle + + Cone slope; // Description of the slope angle + Cone draft; // Description of the draft angle + }; + private: - Double_t fMaxR; + + // -------------------------------------------------- + + Double_t fMaxR; //! + + // ------------- Properties of the lens -------------- + + double fF; // Focal length + double fR; // Radius of lens + double fW; // Width of groove + double fH; // Thickness of lens + double fN; //! Reference refractive index for lens design + double fLambda; // Wavelength [nm] corresponding to fN (at which lens is optimized) + double fVc; //! Velocity of light within the lens material [cm/ns] + double fPSF; // Measure for the surface roughness + + std::vector fGrooves; //! Collection of all grooves + + // ----------- Properties of the material ------------- + + MSpline3 fAbsorptionLength; // Spline storing the absorption length vs wavelength + + // ---------------------------------------------------- + + bool fSlopeAbsorption; //! Absorb rays which hit the slope surface + bool fDraftAbsorption; //! Absorb rays which hit the draft surface + bool fBottomReflection; //! Absorb rays which would be reflected at the exit surface + bool fDisableMultiEntry; //! Absorb upgoing rays inside the material which do not originate from the exit surface + bool fFresnelReflection; //! Disable Fresnel reflection + UInt_t fMinHits; //! Minimum number of surface contacts to define a successfull passage + UInt_t fMaxHits; //! Maximum number of surface contacts to define a successfull passage (0=unlimited) void InitMaxR(); + void InitGeometry(double maxr, double width, double N0, double F, double d); + bool Transmission(double dt, double lambda) const; + + double CalcIntersection(const MQuaternion &p, const MQuaternion &u, const Cone &cone) const; + int FindPeak(size_t i, const MQuaternion &p, const MQuaternion &u) const; + + //int EnterGroove(int surface, double n0, double lambda, MQuaternion &pos, MQuaternion &dir) const; + //int LeavePeak(int surface, double n0, double lambda, MQuaternion &pos, MQuaternion &dir, double T0) const; + + int EnterGroove(int surface, double n0, MQuaternion &pos, MQuaternion &dir) const; + int LeavePeak(int surface, double n0, MQuaternion &pos, MQuaternion &dir, double T0) const; public: @@ -25,14 +110,51 @@ Int_t ExecuteOptics(MQuaternion &p, MQuaternion &u, const Short_t &) const; + Int_t TraceRay(std::vector &vec, MQuaternion &p, MQuaternion &u, const Short_t &wavelength, bool verbose=false) const; - Bool_t IsValid() const { return kTRUE; } + Bool_t IsValid() const { return fGrooves.size(); } + + // ----------------------------------------------------------- + + void SetPSF(const double &psf) { fPSF = psf; } + void DefineLens(double F=50.21, double D=54.92, double w=0.01, double h=0.25, double lambda=546); + + Int_t ReadTransmission(const TString &file, float thickness, bool correction=true); + + // ----------------------------------------------------------- + + void EnableSlopeAbsorption(bool abs=true) { fSlopeAbsorption = abs; } + void EnableDraftAbsorption(bool abs=true) { fDraftAbsorption = abs; } + void DisableBottomReflection(bool ref=true) { fBottomReflection = !ref; } + void DisableMultiEntry(bool mul=true) { fDisableMultiEntry = mul; } + void DisableFresnelReflection(bool ref=true) { fFresnelReflection = !ref; } + + void SetMinHits(UInt_t min) { fMinHits = min; } + void SetMaxHits(UInt_t max) { fMaxHits = max; } + + // ----------------------------------------------------------- + + const std::vector GetGrooves() const { return fGrooves; } + const Groove &GetGroove(int i) const { return fGrooves[i]; } + const Groove &operator[](int i) const { return fGrooves[i]; } + size_t GetNumGrooves() const { return fGrooves.size(); } // ----------------------------------------------------------- static double RefractiveIndex(double lambda); - static double SlopeAngle(double r, double F, double n); + static double SlopeAngle(double r, double F, double n, double d); static double DraftAngle(double r); - ClassDef(MFresnelLens, 1) // Parameter container storing the description of a lens + static double SlopeAngleParabolic(double r, double F, double n0, double n1, double d); + static double SlopeAngleAspherical(double r); + static double SlopeAngleOptimized(double r, double F, double n); + + double SlopeAngle(const double &r) const + { + return SlopeAngle(r, fF, fN, fH); + } + + // ----------------------------------------------------------- + + ClassDef(MFresnelLens, 1) // Parameter container storing the description of a fresnel lens };