sysrap/tests/S4OpBoundaryProcessTest.cc

0001 /**
0002 S4OpBoundaryProcessTest.cc
0003 ============================
0004
0005 Note that the test momentum here is directly
0006 against the normal... changing that leads to
0007 very different distribs of normals.
0008
0009 Added use of S4Random to use precooked randoms
0010 to allow alignment
0011
0012 **/
0013
0014 #include "ssys.h"
0015 #include "spath.h"
0016 #include "NP.hh"
0017 #include "S4OpBoundaryProcess.h"
0018 #include "S4Random.h"
0019
0020 struct Chk
0021 {
0022     Chk();
0023     std::string desc() const ;
0024
0025     const char* FOLD ;
0026     const char* CHECK ;
0027     NP*         a  ;
0028
0029     void run();
0030     void SmearNormal(int chk, double value);
0031
0032 };
0033
0034
0035 Chk::Chk()
0036     :
0037     FOLD(ssys::getenvvar("FOLD")),
0038     CHECK(FOLD ? spath::Basename(FOLD) : nullptr),
0039     a(nullptr)
0040 {
0041 }
0042
0043 std::string Chk::desc() const
0044 {
0045     std::stringstream ss ;
0046     ss << "Chk::desc"
0047        << std::endl
0048        << " FOLD " << ( FOLD ? FOLD : "-" )
0049        << std::endl
0050        << " CHECK " << ( CHECK ? CHECK : "-" )
0051        << std::endl
0052        ;
0053     std::string str = ss.str();
0054     return str ;
0055 }
0056
0057
0058 void Chk::SmearNormal(int chk, double value)
0059 {
0060     S4Random r ;
0061
0062     G4ThreeVector Momentum(0., 0., -1. );
0063     G4ThreeVector Normal(  0., 0.,  1. );
0064
0065     const int N = ssys::getenvint("NUM", 1000) ;
0066     int ni = N ;
0067     int nj = 4 ;
0068
0069     NP* a = NP::Make<double>( ni, nj );
0070     double* aa = a->values<double>();
0071
0072     a->set_meta<double>("value", value );
0073     a->set_meta<std::string>("valuename", chk == 0 ? "sigma_alpha" : "polish" );
0074     a->set_meta<std::string>("source", "S4OpBoundaryProcessTest.sh") ;
0075
0076     for(int i=0 ; i < ni ; i++)
0077     {
0078         r.setSequenceIndex(i);   // use precooked random streams : so can align
0079
0080         G4ThreeVector FacetNormal( 0., 0., 0.) ;
0081         switch(chk)
0082         {
0083             case 0: FacetNormal = S4OpBoundaryProcess::SmearNormal_SigmaAlpha( Momentum, Normal, value );  break ;
0084             case 1: FacetNormal = S4OpBoundaryProcess::SmearNormal_Polish(     Momentum, Normal, value );  break ;
0085         }
0086
0087         aa[i*4+0] = FacetNormal.x() ;
0088         aa[i*4+1] = FacetNormal.y() ;
0089         aa[i*4+2] = FacetNormal.z() ;
0090         aa[i*4+3] = 0.  ;
0091
0092         if(i < 4)  std::cout
0093             << " i " << i
0094             << " FacetNormal = np.array(["
0095             << std::fixed << std::setw(10) << std::setprecision(5) << FacetNormal.x() << ","
0096             << std::fixed << std::setw(10) << std::setprecision(5) << FacetNormal.y() << ","
0097             << std::fixed << std::setw(10) << std::setprecision(5) << FacetNormal.z() << "])"
0098             << std::endl
0099             ;
0100     }
0101     a->save("$FOLD/q.npy" );
0102 }
0103
0104 void Chk::run()
0105 {
0106     std::cout << desc() << std::endl ;
0107
0108     if(     strcmp(CHECK, "smear_normal_sigma_alpha")==0) SmearNormal(0, 0.1) ;
0109     else if(strcmp(CHECK, "smear_normal_polish")==0)      SmearNormal(1, 0.8) ;
0110     else
0111     {
0112         std::cerr
0113             << "Chk::run"
0114             << " CHECK " << ( CHECK ? CHECK : "-" ) << " IS UNHANDLED "
0115             << std::endl
0116             ;
0117     }
0118 }
0119
0120
0121 int main()
0122 {
0123     Chk chk ;
0124     chk.run();
0125
0126     return 0 ;
0127 }