eic-opticks/sysrap/SThetaCut.cc

0001 #include <csignal>
0002 #include "SLOG.hh"
0003 #include "scuda.h"
0004 #include "squad.h"
0005 #include "SThetaCut.hh"
0006
0007 const plog::Severity SThetaCut::LEVEL = SLOG::EnvLevel("SThetaCut", "DEBUG" );
0008
0009 /**
0010 SThetaCut::PrepareParam
0011 ------------------------
0012
0013 This is used by
0014
0015 * npy/NThetaCut.cpp
0016 * CSG/CSGNode.cc
0017
0018
0019 The input parameters: startTheta, deltaTheta configure two polar angles in range [0., 1.] (units of pi)
0020
0021 * startTheta
0022 * startTheta+deltaTheta
0023
0024 A do nothing thetacut "cutter" would have have startTheta=0. deltaTheta=1.
0025
0026
0027            Z
0028
0029            0.0
0030        \   :   /
0031         \  :  /
0032          \ : /
0033           \:/
0034            O - - -  0.5    X
0035           /:\
0036          / : \
0037         /  :  \
0038        /   :   \
0039            1.0
0040
0041
0042 +---------------------------+----------------------+---------------------------------+-------------------------------+
0043 | (startTheta, deltaTheta)  | (theta0, theta1)     |  Notes                          |  Cutaway                      |
0044 +===========================+======================+=================================+===============================+
0045 |  (0.0, 1.0)               |  (0.0, 1.0)          |  Full range of theta            |                               |
0046 +---------------------------+----------------------+---------------------------------+-------------------------------+
0047 |  (0.1, 0.9)               |  (0.1, 1.0)          |  top hole around zenith         |  0.0->0.1                     |
0048 +---------------------------+----------------------+---------------------------------+-------------------------------+
0049 |  (0.0, 0.9)               |  (0.0, 0.9)          |  bot hole around nadir          |  0.9->1.0                     |
0050 +---------------------------+----------------------+---------------------------------+-------------------------------+
0051 |  (0.25, 0.50)             |  (0.25, 0.75)        |  butterfly or bowtie            |  0.0->0.25, 0.75->1.0         |
0052 +---------------------------+----------------------+---------------------------------+-------------------------------+
0053 |  (0.45, 0.10)             |  (0.45, 0.55)        |  twisted belt                   |  0.0->0.45, 0.55->1.0         |
0054 +---------------------------+----------------------+---------------------------------+-------------------------------+
0055
0056 **/
0057
0058
0059 void SThetaCut::PrepareParam( quad& q0, quad& q1, double startTheta_pi, double deltaTheta_pi )
0060 {
0061     double theta0_pi = startTheta_pi  ;
0062     double theta1_pi = startTheta_pi + deltaTheta_pi ;
0063     bool has_thetacut = theta0_pi > 0. || theta1_pi < 1. ;
0064     assert(has_thetacut);
0065     if(!has_thetacut) std::raise(SIGINT);
0066
0067     assert( theta0_pi >= 0. && theta0_pi <= 1.) ;
0068     assert( theta1_pi >= 0. && theta1_pi <= 1.) ;
0069
0070     bool plane_theta0 = theta0_pi == 0.5 ;
0071     bool plane_theta1 = theta1_pi == 0.5 ;
0072     // one of the cones can degenerate to a plane for angle 0.5
0073     // if they both do this then then nothing is left
0074
0075     const double pi = M_PIf ;
0076     double theta0 = theta0_pi*pi ;
0077     double theta1 = theta1_pi*pi ;
0078
0079     double cosTheta0 = std::cos(theta0);
0080     double sinTheta0 = std::sin(theta0);
0081     double tanTheta0 = std::tan(theta0);
0082
0083     double cosTheta1 = std::cos(theta1);
0084     double sinTheta1 = std::sin(theta1);
0085     double tanTheta1 = std::tan(theta1);
0086
0087     //double cosTheta0Sign = cosTheta0/std::abs(cosTheta0) ;
0088     //double cosTheta1Sign = cosTheta1/std::abs(cosTheta1) ;
0089
0090     double tan2Theta0 = tanTheta0*tanTheta0 ;
0091     double tan2Theta1 = tanTheta1*tanTheta1 ;
0092
0093
0094     q0.f.x = float(plane_theta0 ? 0.  : cosTheta0) ;
0095     q0.f.y = float(plane_theta0 ? 1.0 : sinTheta0) ;
0096     q0.f.z = float(plane_theta1 ? 0.  : cosTheta1) ;
0097     q0.f.w = float(plane_theta1 ? 1.0 : sinTheta1) ;
0098
0099     q1.f.x = float(plane_theta0 ? 0. : tan2Theta0 );  // skipping the infinity
0100     q1.f.y = float(plane_theta1 ? 0. : tan2Theta1 );
0101
0102
0103     LOG(LEVEL)
0104         << " startTheta_pi " << startTheta_pi
0105         << " deltaTheta_pi " << deltaTheta_pi
0106         << " plane_theta0 " << plane_theta0
0107         << " plane_theta1 " << plane_theta1
0108         << " theta0 " << theta0
0109         << " theta1 " << theta1
0110         ;
0111
0112     LOG(LEVEL)
0113         << " cosTheta0 " << cosTheta0
0114         << " sinTheta0 " << sinTheta0
0115         << " cosTheta1 " << cosTheta1
0116         << " sinTheta1 " << sinTheta1
0117         << " tan2Theta0 " << tan2Theta0
0118         << " tan2Theta1 " << tan2Theta1
0119         ;
0120
0121
0122 }
0123