eic-opticks/sysrap/AngleCut.cc

0001 #include <cassert>
0002 #include <cmath>
0003 #include <csignal>
0004
0005 #include "AngleCut.hh"
0006 #include "SLOG.hh"
0007 #include "SMath.hh"
0008 #include "scuda.h"
0009 #include "squad.h"
0010
0011 const plog::Severity SPhiCut::LEVEL = SLOG::EnvLevel("SPhiCut", "DEBUG");
0012
0013 /**
0014 SPhiCut::PrepareParam
0015 ------------------------
0016
0017 This is used by
0018
0019 * npy/NPhiCut.cpp
0020 * CSG/CSGNode.cc
0021
0022
0023 The input parameters: startPhi, deltaPhi configure two azimuthal angles in range [0., 2.] (units of pi)
0024
0025 * startPhi
0026 * startPhi+deltaPhi
0027
0028 A do nothing phicut shape would have have startPhi=0. deltaPhi=2.
0029
0030
0031                Y
0032
0033               0.5
0034                :   /
0035                :  /
0036                : /
0037                :/
0038      1.0- - -  O - - -  0.0    X
0039                :\
0040                : \
0041                :  \
0042                :   \
0043
0044               1.5
0045
0046
0047 +---------------------------+----------------------+---------------------------------+-------------------------------+
0048 | (startPhi, deltaPhi)      | (phi0, phi1)         |  Notes                          |  Cutaway                      |
0049 +===========================+======================+=================================+===============================+
0050 |  (0.0, 2.0)               |  (0.0, 2.0)          |  Full range of phi              |                               |
0051 +---------------------------+----------------------+---------------------------------+-------------------------------+
0052
0053 **/
0054
0055 void SPhiCut::PrepareParam(quad& q0, double startPhi_pi, double deltaPhi_pi)
0056 {
0057     double phi0_pi = startPhi_pi;
0058     double phi1_pi = startPhi_pi + deltaPhi_pi;
0059
0060     double cosPhi0 = SMath::cos_pi(phi0_pi);
0061     double sinPhi0 = SMath::sin_pi(phi0_pi);
0062     double cosPhi1 = SMath::cos_pi(phi1_pi);
0063     double sinPhi1 = SMath::sin_pi(phi1_pi);
0064
0065     q0.f.x = float(cosPhi0);
0066     q0.f.y = float(sinPhi0);
0067     q0.f.z = float(cosPhi1);
0068     q0.f.w = float(sinPhi1);
0069
0070     LOG(LEVEL)
0071         << " startPhi_pi " << startPhi_pi
0072         << " deltaPhi_pi " << deltaPhi_pi
0073         << " phi0_pi " << phi0_pi
0074         << " phi1_pi " << phi1_pi
0075         << " cosPhi0 " << cosPhi0
0076         << " sinPhi0 " << sinPhi0
0077         << " cosPhi1 " << cosPhi1
0078         << " sinPhi1 " << sinPhi1;
0079 }
0080
0081 const plog::Severity SThetaCut::LEVEL = SLOG::EnvLevel("SThetaCut", "DEBUG");
0082
0083 /**
0084 SThetaCut::PrepareParam
0085 ------------------------
0086
0087 This is used by
0088
0089 * npy/NThetaCut.cpp
0090 * CSG/CSGNode.cc
0091
0092
0093 The input parameters: startTheta, deltaTheta configure two polar angles in range [0., 1.] (units of pi)
0094
0095 * startTheta
0096 * startTheta+deltaTheta
0097
0098 A do nothing thetacut "cutter" would have have startTheta=0. deltaTheta=1.
0099
0100
0101            Z
0102
0103            0.0
0104        \   :   /
0105         \  :  /
0106          \ : /
0107           \:/
0108            O - - -  0.5    X
0109           /:\
0110          / : \
0111         /  :  \
0112        /   :   \
0113            1.0
0114
0115
0116 +---------------------------+----------------------+---------------------------------+-------------------------------+
0117 | (startTheta, deltaTheta)  | (theta0, theta1)     |  Notes                          |  Cutaway                      |
0118 +===========================+======================+=================================+===============================+
0119 |  (0.0, 1.0)               |  (0.0, 1.0)          |  Full range of theta            |                               |
0120 +---------------------------+----------------------+---------------------------------+-------------------------------+
0121 |  (0.1, 0.9)               |  (0.1, 1.0)          |  top hole around zenith         |  0.0->0.1                     |
0122 +---------------------------+----------------------+---------------------------------+-------------------------------+
0123 |  (0.0, 0.9)               |  (0.0, 0.9)          |  bot hole around nadir          |  0.9->1.0                     |
0124 +---------------------------+----------------------+---------------------------------+-------------------------------+
0125 |  (0.25, 0.50)             |  (0.25, 0.75)        |  butterfly or bowtie            |  0.0->0.25, 0.75->1.0         |
0126 +---------------------------+----------------------+---------------------------------+-------------------------------+
0127 |  (0.45, 0.10)             |  (0.45, 0.55)        |  twisted belt                   |  0.0->0.45, 0.55->1.0         |
0128 +---------------------------+----------------------+---------------------------------+-------------------------------+
0129
0130 **/
0131
0132 void SThetaCut::PrepareParam(quad& q0, quad& q1, double startTheta_pi, double deltaTheta_pi)
0133 {
0134     double theta0_pi = startTheta_pi;
0135     double theta1_pi = startTheta_pi + deltaTheta_pi;
0136     bool   has_thetacut = theta0_pi > 0. || theta1_pi < 1.;
0137     assert(has_thetacut);
0138     if (!has_thetacut)
0139         std::raise(SIGINT);
0140
0141     assert(theta0_pi >= 0. && theta0_pi <= 1.);
0142     assert(theta1_pi >= 0. && theta1_pi <= 1.);
0143
0144     bool plane_theta0 = theta0_pi == 0.5;
0145     bool plane_theta1 = theta1_pi == 0.5;
0146     // one of the cones can degenerate to a plane for angle 0.5
0147     // if they both do this then then nothing is left
0148
0149     const double pi = M_PIf;
0150     double       theta0 = theta0_pi * pi;
0151     double       theta1 = theta1_pi * pi;
0152
0153     double cosTheta0 = std::cos(theta0);
0154     double sinTheta0 = std::sin(theta0);
0155     double tanTheta0 = std::tan(theta0);
0156
0157     double cosTheta1 = std::cos(theta1);
0158     double sinTheta1 = std::sin(theta1);
0159     double tanTheta1 = std::tan(theta1);
0160
0161     // double cosTheta0Sign = cosTheta0/std::abs(cosTheta0) ;
0162     // double cosTheta1Sign = cosTheta1/std::abs(cosTheta1) ;
0163
0164     double tan2Theta0 = tanTheta0 * tanTheta0;
0165     double tan2Theta1 = tanTheta1 * tanTheta1;
0166
0167     q0.f.x = float(plane_theta0 ? 0. : cosTheta0);
0168     q0.f.y = float(plane_theta0 ? 1.0 : sinTheta0);
0169     q0.f.z = float(plane_theta1 ? 0. : cosTheta1);
0170     q0.f.w = float(plane_theta1 ? 1.0 : sinTheta1);
0171
0172     q1.f.x = float(plane_theta0 ? 0. : tan2Theta0); // skipping the infinity
0173     q1.f.y = float(plane_theta1 ? 0. : tan2Theta1);
0174
0175     LOG(LEVEL)
0176         << " startTheta_pi " << startTheta_pi
0177         << " deltaTheta_pi " << deltaTheta_pi
0178         << " plane_theta0 " << plane_theta0
0179         << " plane_theta1 " << plane_theta1
0180         << " theta0 " << theta0
0181         << " theta1 " << theta1;
0182
0183     LOG(LEVEL)
0184         << " cosTheta0 " << cosTheta0
0185         << " sinTheta0 " << sinTheta0
0186         << " cosTheta1 " << cosTheta1
0187         << " sinTheta1 " << sinTheta1
0188         << " tan2Theta0 " << tan2Theta0
0189         << " tan2Theta1 " << tan2Theta1;
0190 }