eic-opticks/ana/twhite.py

0001 #!/usr/bin/env python
0002 #
0003 # Copyright (c) 2019 Opticks Team. All Rights Reserved.
0004 #
0005 # This file is part of Opticks
0006 # (see https://bitbucket.org/simoncblyth/opticks).
0007 #
0008 # Licensed under the Apache License, Version 2.0 (the "License");
0009 # you may not use this file except in compliance with the License.
0010 # You may obtain a copy of the License at
0011 #
0012 #   http://www.apache.org/licenses/LICENSE-2.0
0013 #
0014 # Unless required by applicable law or agreed to in writing, software
0015 # distributed under the License is distributed on an "AS IS" BASIS,
0016 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
0017 # See the License for the specific language governing permissions and
0018 # limitations under the License.
0019 #
0020
0021 """
0022 twhite.py: Wavelength Distribution Check
0023 ============================================
0024
0025 Creates plot comparing simulated photon wavelength spectrum
0026 from :doc:`../tests/twhite` against blackbody expectation.
0027
0028 This is checking the *source_lookup* implementation and
0029 the inverse CDF *source_texture* that it uses.
0030
0031 .. code-block:: cpp
0032
0033     ## optixrap/cu/wavelength_lookup.h
0034
0035     014 rtTextureSampler<float, 2>  source_texture ;
0036      15 rtDeclareVariable(float4, source_domain, , );
0037      ..
0038      41 static __device__ __inline__ float source_lookup(float u)
0039      42 {
0040      43     float ui = u/source_domain.z + 0.5f ;
0041      44     return tex2D(source_texture, ui, 0.5f );  // line 0
0042      45 }
0043      46
0044      47 static __device__ __inline__ void source_check()
0045      48 {
0046      49     float nm_a = source_lookup(0.0f);
0047      50     float nm_b = source_lookup(0.5f);
0048      51     float nm_c = source_lookup(1.0f);
0049      52     rtPrintf("source_check nm_a %10.3f %10.3f %10.3f  \n",  nm_a, nm_b, nm_c );
0050      53 }
0051
0052     ## optixrap/cu/torchstep.h
0053
0054     241 __device__ void
0055     242 generate_torch_photon(Photon& p, TorchStep& ts, RNG &rng)
0056     243 {
0057     244       p.wavelength = ts.wavelength > 50.f ? ts.wavelength : source_lookup(curand_uniform(&rng));  // Planck black body source 6500K standard illuminant
0058     245
0059
0060
0061 See Also
0062 ----------
0063
0064 * :doc:`source_debug`
0065
0066
0067 """
0068
0069 import os, sys, logging, numpy as np
0070 log = logging.getLogger(__name__)
0071
0072 import matplotlib.pyplot as plt
0073 from mpl_toolkits.mplot3d import Axes3D
0074
0075 from opticks.ana.base import opticks_main
0076 from opticks.ana.evt import Evt
0077 from opticks.ana.planck import planck
0078
0079
0080
0081 if __name__ == '__main__':
0082     plt.ion()
0083     args = opticks_main(tag="1", det="white", src="torch")
0084
0085     ## tag = "1"   ## dont have any tag 1 anymore
0086     ## tag = "15"     ## so added tag 15,16 to ggv-rainbow with wavelength=0 which is default black body
0087
0088     try:
0089         evt = Evt(tag=args.tag, det=args.det, src=args.src, args=args )
0090     except IOError as err:
0091         log.fatal(err)
0092         sys.exit(args.mrc)
0093
0094
0095     if not evt.valid:
0096        log.fatal("failed to load evt %s " % repr(args))
0097        sys.exit(1)
0098
0099
0100     wl = evt.wl
0101     w0 = evt.recwavelength(0)
0102
0103     w = wl
0104     #w = w0
0105
0106     wd = np.linspace(60,820,256) - 1.
0107     # reduce bin edges by 1nm to avoid aliasing artifact in the histogram
0108
0109     mid = (wd[:-1]+wd[1:])/2.     # bin middle
0110
0111     pl = planck(mid, 6500.)
0112     pl /= pl.sum()
0113
0114     counts, edges = np.histogram(w, bins=wd )
0115     fcounts = counts.astype(np.float32)
0116     fcounts  /= fcounts.sum()
0117
0118
0119     plt.close()
0120
0121     plt.plot( edges[:-1], fcounts, drawstyle="steps-mid")
0122
0123     plt.plot( mid,  pl )
0124
0125     plt.axis( [w.min() - 100, w.max() + 100, 0, fcounts.max()*1.1 ])
0126
0127     #plt.hist(w, bins=256)   # 256 is number of unique wavelengths (from record compression)
0128
0129