eic-opticks/ana/rayleigh.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 rayleigh.py
0023 =============================================
0024
0025 Without selection scatter distrib plots from
0026 arrays created by:
0027
0028 * optixrap/tests/ORayleighTest.cc
0029 * cfg4/tests/OpRayleighTest.cc
0030
0031
0032 """
0033 import os, sys, logging, numpy as np
0034 log = logging.getLogger(__name__)
0035
0036 import matplotlib.pyplot as plt
0037 from opticks.ana.nbase import vnorm, costheta_
0038
0039
0040 OLDMOM,OLDPOL,NEWMOM,NEWPOL = 0,1,2,3
0041 X,Y,Z=0,1,2
0042
0043
0044
0045 def dotmom_(a):
0046     oldmom = a[:,OLDMOM,:3]
0047     newmom = a[:,NEWMOM,:3]
0048     dotmom = costheta_(oldmom,newmom)
0049     return dotmom
0050
0051 def dotpol_(a):
0052     oldpol = a[:,OLDPOL,:3]
0053     newpol = a[:,NEWPOL,:3]
0054     dotpol = costheta_(oldpol,newpol)
0055     return dotpol
0056
0057
0058 if __name__ == '__main__':
0059
0060     aa = np.load(os.path.expandvars("$TMP/RayleighTest/ok.npy"))
0061     bb = np.load(os.path.expandvars("$TMP/RayleighTest/cfg4.npy"))
0062
0063     bins = 100
0064     nx = 4
0065     ny = 2
0066
0067     qwns = [
0068          (1,aa[:,NEWMOM,X],bb[:,NEWMOM,X],"momx"),
0069          (2,aa[:,NEWMOM,Y],bb[:,NEWMOM,Y],"momy"),
0070          (3,aa[:,NEWMOM,Z],bb[:,NEWMOM,Z],"momz"),
0071          (4,dotmom_(aa)   ,dotmom_(bb)   ,"dotmom"),
0072
0073          (5,aa[:,NEWPOL,X],bb[:,NEWPOL,X],"polx"),
0074          (6,aa[:,NEWPOL,Y],bb[:,NEWPOL,Y],"poly"),
0075          (7,aa[:,NEWPOL,Z],bb[:,NEWPOL,Z],"polz"),
0076          (8,dotpol_(aa)   ,dotpol_(bb)   ,"dotpol"),
0077            ]
0078
0079     for i,a,b,label in qwns:
0080         plt.subplot(ny, nx, i)
0081         plt.hist(a, bins=bins, histtype="step", label=label)
0082         plt.hist(b, bins=bins, histtype="step", label=label)
0083     pass
0084     plt.show()
0085
0086