u4/tests/G4OpBoundaryProcess_GetFacetNormal_Test.cc

0001 // ./G4OpBoundaryProcess_GetFacetNormal_Test.sh
0002
0003 #include <iostream>
0004 #include "G4ThreeVector.hh"
0005 #include "Randomize.hh"
0006 #include "NPX.h"
0007
0008 /**
0009 GetFacetNormal
0010 ----------------
0011
0012 Only Momentum_FacetNormal < 0. escapes the while loop
0013 ie FacetNormal must end up against the Momentum.
0014
0015 Presumably this is assuming the initial Momentum is against the Normal.
0016 It would make more sense to constrain FacetNormal to
0017 staying with the same orientation as the initial one
0018 not just assuming a particular orientation.
0019
0020
0021 Smearing with polish 0.999 results in something effectively indistingishable from specular,
0022 so it just consumes randoms and does little else.
0023
0024   : t.FacetNormal                                      :            (1000, 3) : 0:03:32.992163
0025
0026  min_stamp : 2023-03-03 17:52:04.431359
0027  max_stamp : 2023-03-03 17:52:04.431359
0028  dif_stamp : 0:00:00
0029  age_stamp : 0:03:32.992163
0030 [[-0.000008198991 -0.000774490262  0.999999700049]
0031  [-0.000053926407  0.000675428761  0.999999770444]
0032  [-0.000123759488 -0.000047115051  0.999999991232]
0033  [-0.000732900199 -0.000029525774  0.999999730993]
0034  [-0.000347653389 -0.000362541626  0.99999987385 ]
0035  [-0.000390777914  0.000540491843  0.999999777581]
0036  [-0.000222860429 -0.000332362899  0.999999919934]
0037  [-0.000550520639 -0.000345155409  0.999999788897]
0038  [-0.000083752166 -0.000284001987  0.999999956164]
0039  [-0.000172539198 -0.000299258882  0.999999940337]
0040  [-0.000378540935  0.000299374342  0.999999883541]
0041  [-0.000554646197  0.000462350078  0.9999997393  ]
0042  [-0.000523868695 -0.000421828994  0.999999773811]
0043  [-0.00046533652   0.0006765364    0.99999966288 ]
0044  [-0.000070471159  0.000012100613  0.999999997444]
0045  [-0.000041047925  0.0000830588    0.999999995708]
0046  [-0.000027495856  0.000620312374  0.999999807228]
0047  [-0.000708336559 -0.000352951962  0.999999686842]
0048  [-0.00005016816   0.000631667703  0.99999979924 ]
0049  [-0.000074296972  0.000283778564  0.999999956975]
0050
0051 **/
0052
0053
0054 G4ThreeVector GetFacetNormal(G4ThreeVector& smear, const G4ThreeVector& Momentum, const G4ThreeVector& Normal, G4double polish )
0055 {
0056     G4ThreeVector FacetNormal;
0057     G4double Momentum_FacetNormal = 0. ;
0058
0059     if (polish < 1.0)
0060     {
0061         do {
0062             do
0063             {
0064                 smear.setX(2.*G4UniformRand()-1.0);
0065                 smear.setY(2.*G4UniformRand()-1.0);
0066                 smear.setZ(2.*G4UniformRand()-1.0);
0067             }
0068             while (smear.mag()>1.0);     // random vector within unit sphere (not normalized)
0069
0070             smear = (1.-polish) * smear;   // scale it down (greatly for polish 0.999)
0071
0072             FacetNormal = Normal + smear;  // perturb the Normal by the smear
0073
0074             Momentum_FacetNormal = Momentum * FacetNormal ;
0075
0076         }
0077         while (Momentum_FacetNormal >= 0.0);
0078
0079         FacetNormal = FacetNormal.unit();
0080     }
0081     else
0082     {
0083         FacetNormal = Normal;
0084     }
0085     return FacetNormal ;
0086 }
0087
0088
0089 struct FacetNormal_Smear
0090 {
0091     G4ThreeVector FacetNormal ;
0092     G4ThreeVector Smear ;
0093 };
0094
0095 int main()
0096 {
0097     G4double      polish = 0.999 ;
0098
0099     G4ThreeVector Momentum(1,-1,0);
0100     G4ThreeVector Normal(0,0,1) ;
0101     G4ThreeVector Meta(polish, 0, 0);
0102
0103     Momentum = Momentum.unit();
0104     Normal = Normal.unit() ;
0105
0106     std::vector<G4ThreeVector> mm(3) ;
0107     mm[0] = Momentum ;
0108     mm[1] = Normal ;
0109     mm[2] = Meta ;
0110
0111     const int N=1000 ;
0112     std::vector<FacetNormal_Smear> nn(N) ;
0113
0114     for(int i=0 ; i < int(nn.size()) ; i++ )
0115     {
0116         FacetNormal_Smear& fns = nn[i] ;
0117         fns.FacetNormal = GetFacetNormal(fns.Smear, Momentum, Normal, polish ) ;
0118     }
0119
0120     NP* a = NPX::ArrayFromVec<double,FacetNormal_Smear>(nn,2,3) ;
0121     a->save("$FOLD/FacetNormal.npy" );
0122
0123     NP* b = NPX::ArrayFromVec<double,G4ThreeVector>(mm) ;
0124     b->save("$FOLD/Meta.npy" );
0125
0126     return 0 ;
0127 }
0128
0129