sysrap/tests/S4OpBoundaryProcess.h

0001 #pragma once
0002
0003 #include "G4ThreeVector.hh"
0004
0005 struct S4OpBoundaryProcess
0006 {
0007     static G4ThreeVector SmearNormal_SigmaAlpha(const G4ThreeVector& Momentum, const G4ThreeVector&  Normal, G4double sigma_alpha );
0008     static G4ThreeVector SmearNormal_Polish(    const G4ThreeVector& Momentum, const G4ThreeVector&  Normal, G4double polish      );
0009 };
0010
0011
0012 #include "Randomize.hh"
0013 #include <CLHEP/Units/PhysicalConstants.h>
0014
0015 /**
0016 S4OpBoundaryProcess::SmearNormal_SigmaAlpha
0017 --------------------------------------------
0018
0019 Extract from G4OpBoundaryProcess::GetFacetNormal to investigate how to
0020 do the equivalent with CUDA.
0021
0022 HMM: can implement all this in CUDA : the difficulty is more with
0023 how to implement switching to alternative "ground" rather than "polished"
0024 handling. Maybe best way to switch by expanding the ems enum used in qsim::propagate::
0025
0026     2025     if( command == BOUNDARY )
0027     2026     {
0028     2027         const int& ems = ctx.s.optical.y ;
0029     ....
0030     2040         if( ems == smatsur_NoSurface )
0031     2041         {
0032     2042             command = propagate_at_boundary( flag, rng, ctx ) ;
0033     2043         }
0034     2044         else if( ems == smatsur_Surface )
0035     2045         {
0036     2046             command = propagate_at_surface( flag, rng, ctx ) ;
0037     2047         }
0038
0039 Or branch based on ems or some other surface type enum within::
0040
0041     qsim::propagate_at_surface
0042
0043 * HMM: probably best to defer this until find a corresponding discrepancy ?
0044   So can change just the thing to get match
0045
0046
0047 Q1:Where to get sigma_alpha and polish from ?
0048 A1:U4Surface::MakeFold plants ModelValue, etc.. into surface metadata::
0049
0050     288         G4double ModelValue = theModel == glisur ? os->GetPolish() : os->GetSigmaAlpha() ;
0051     289         assert( ModelValue >= 0. && ModelValue <= 1. );
0052     ...
0053     302         sub->set_meta<int>("Type", theType) ;
0054     303         sub->set_meta<int>("Model", theModel) ;
0055     304         sub->set_meta<int>("Finish", theFinish) ;
0056     305         sub->set_meta<double>("ModelValue", ModelValue ) ;
0057
0058 The value percentage is in optical buffer .w already
0059 could trivially change to setting float32 value into
0060 the mostly int optical buffer.::
0061
0062     In [8]: f.bnd_names[4]
0063     Out[8]: 'Water/NNVTMaskOpticalSurface//CDReflectorSteel'
0064
0065     In [9]: f.optical[4]
0066     Out[9]:
0067     array([[ 3,  0,  0,  0],
0068            [27,  2,  3, 20],
0069            [ 0,  1,  0,  0],
0070            [ 2,  0,  0,  0]], dtype=int32)
0071
0072
0073 Q2: How to switch on FacetNormal smearing with SmearNormal_SigmaAlpha/SmearNormal_Polish ?
0074     How to pick between them ?
0075
0076
0077
0078 **/
0079
0080 inline G4ThreeVector S4OpBoundaryProcess::SmearNormal_SigmaAlpha(const G4ThreeVector& Momentum, const G4ThreeVector&  Normal, G4double sigma_alpha )
0081 {
0082     if (sigma_alpha == 0.0) return Normal;
0083     G4ThreeVector FacetNormal ;
0084
0085     G4double alpha;
0086     G4double f_max = std::min(1.0,4.*sigma_alpha);
0087
0088     G4double phi, SinAlpha, CosAlpha, SinPhi, CosPhi, unit_x, unit_y, unit_z;
0089     G4ThreeVector tmpNormal;
0090
0091     do {
0092        do {
0093           alpha = G4RandGauss::shoot(0.0,sigma_alpha);
0094           // Loop checking, 13-Aug-2015, Peter Gumplinger
0095        } while (G4UniformRand()*f_max > std::sin(alpha) || alpha >= CLHEP::halfpi );
0096
0097        phi = G4UniformRand()*CLHEP::twopi;
0098
0099        SinAlpha = std::sin(alpha);
0100        CosAlpha = std::cos(alpha);
0101        SinPhi = std::sin(phi);
0102        CosPhi = std::cos(phi);
0103
0104        unit_x = SinAlpha * CosPhi;
0105        unit_y = SinAlpha * SinPhi;
0106        unit_z = CosAlpha;
0107
0108        FacetNormal.setX(unit_x);
0109        FacetNormal.setY(unit_y);
0110        FacetNormal.setZ(unit_z);
0111
0112        tmpNormal = Normal;
0113        // HUH: argument to rotateUz and Normal are both const refs
0114        // so tmpNormal is entirely pointless ?
0115
0116        FacetNormal.rotateUz(tmpNormal);
0117        // Loop checking, 13-Aug-2015, Peter Gumplinger
0118     } while (Momentum * FacetNormal >= 0.0);
0119     return FacetNormal ;
0120 }
0121
0122 inline G4ThreeVector S4OpBoundaryProcess::SmearNormal_Polish(     const G4ThreeVector& Momentum, const G4ThreeVector&  Normal, G4double polish )
0123 {
0124     if (polish == 1.0) return Normal;
0125     G4ThreeVector FacetNormal ;
0126
0127     do {
0128        G4ThreeVector smear;
0129        do {
0130           smear.setX(2.*G4UniformRand()-1.0);
0131           smear.setY(2.*G4UniformRand()-1.0);
0132           smear.setZ(2.*G4UniformRand()-1.0);
0133           // Loop checking, 13-Aug-2015, Peter Gumplinger
0134        } while (smear.mag()>1.0);
0135        smear = (1.-polish) * smear;
0136        FacetNormal = Normal + smear;
0137        // Loop checking, 13-Aug-2015, Peter Gumplinger
0138     } while (Momentum * FacetNormal >= 0.0);
0139     FacetNormal = FacetNormal.unit();
0140     return FacetNormal ;
0141 }
0142