eic-opticks/ana/ckn.py

0001 #!/usr/bin/env python
0002 """
0003 ckn.py : reproduce G4Cerenkov_modified::GetAverageNumberOfPhotons
0004 ====================================================================
0005
0006 ::
0007
0008     ipython -i ckn.py
0009
0010 ::
0011
0012     cp /tmp/blyth/opticks/ana/ckn/* /Users/blyth/simoncblyth.bitbucket.io/env/presentation/ana/ckn/
0013
0014
0015
0016 """
0017 import os, logging, numpy as np
0018 import matplotlib.pyplot as plt
0019 from opticks.ana.main import opticks_main
0020 from opticks.ana.key import keydir
0021 from opticks.ana.nload import np_load
0022
0023 #from scipy import integrate
0024
0025
0026 log = logging.getLogger(__name__)
0027
0028 class CKN(object):
0029     """
0030     Reproduces the G4Cerenkov Frank-Tamm integration to give average number of photons
0031     for a BetaInverse and RINDEX profile.
0032     """
0033     FOLD = os.path.expandvars("/tmp/$USER/opticks/ana/ckn")
0034     kd = keydir()
0035     rindex_path = os.path.join(kd, "GScintillatorLib/LS_ori/RINDEX.npy")
0036
0037     def __init__(self):
0038         ri = np.load(self.rindex_path)
0039         ri[:,0] *= 1e6    # into eV
0040         ri_ = lambda ev:np.interp( ev, ri[:,0], ri[:,1] )
0041
0042         self.ri = ri
0043         self.ri_ = ri_
0044         self.BuildThePhysicsTable()
0045         self.BuildThePhysicsTable_2()
0046         assert np.allclose( self.cai, self.cai2 )
0047         self.paths = []
0048
0049     def BuildThePhysicsTable(self, dump=False):
0050         """
0051         See G4Cerenkov_modified::BuildThePhysicsTable
0052
0053
0054         This is applying the composite trapezoidal rule to do a
0055         numerical energy integral  of  n^(-2) = 1./(ri[:,1]*ri[:,1])
0056         """
0057         ri = self.ri
0058         en = ri[:,0]
0059
0060         ir2 = 1./(ri[:,1]*ri[:,1])
0061         mir2 = 0.5*(ir2[1:] + ir2[:-1])
0062
0063         de = en[1:] - en[:-1]
0064
0065         assert len(mir2) == len(ri) - 1       # averaging points looses one value
0066         mir2_de = mir2*de
0067
0068         cai = np.zeros(len(ri))               # leading zero regains one value
0069         np.cumsum(mir2_de, out=cai[1:])
0070
0071         if dump:
0072             print("cai", cai)
0073         pass
0074
0075         self.cai = cai
0076         self.ir2 = ir2
0077         self.mir2 = mir2
0078         self.de = de
0079         self.mir2_de = mir2_de
0080
0081
0082     def BuildThePhysicsTable_2(self, dump=False):
0083         """
0084         np.trapz does the same thing as above : applying composite trapezoidal integration
0085
0086         https://numpy.org/doc/stable/reference/generated/numpy.trapz.html
0087         """
0088         ri = self.ri
0089         en = ri[:,0]
0090         ir2 = 1./(ri[:,1]*ri[:,1])
0091
0092         cai2 = np.zeros(len(ri))
0093         for i in range(len(ri)):
0094             cai2[i] = np.trapz( ir2[:i+1], en[:i+1] )
0095         pass
0096         self.cai2 = cai2
0097
0098         if dump:
0099             print("cai2", cai2)
0100         pass
0101
0102
0103     @classmethod
0104     def BuildThePhysicsTable_s2i(cls, ri, BetaInverse, dump=False):
0105         """
0106         No need for a physics table when do integral directly on s2
0107         """
0108         en = ri[:,0]
0109         s2i = np.zeros(len(ri))
0110         for i in range(len(ri)):
0111             s2i[i] = np.trapz( s2[:i+1], en[:i+1] )
0112         pass
0113         return s2i
0114
0115     def GetAverageNumberOfPhotons_s2(self, BetaInverse, charge=1, dump=False, s2cross=False ):
0116         """
0117         Simplfied Alternative to _s2messy following C++ implementation.
0118         Allowed regions are identified by s2 being positive avoiding the need for
0119         separately getting crossings. Instead get the crossings and do the trapezoidal
0120         numerical integration in one pass, improving simplicity and accuracy.
0121
0122         See opticks/examples/Geant4/CerenkovStandalone/G4Cerenkov_modified.cc
0123         """
0124         s2integral = 0.
0125         for i in range(len(self.ri)-1):
0126             en = np.array([self.ri[i,0], self.ri[i+1,0] ])
0127             ri = np.array([self.ri[i,1], self.ri[i+1,1] ])
0128             ct = BetaInverse/ri
0129             s2 = (1.-ct)*(1.+ct)
0130
0131             en_0, en_1 = en
0132             ri_0, ri_1 = ri
0133             s2_0, s2_1 = s2
0134
0135             if s2_0*s2_1 < 0.:
0136                 en_cross_A = (s2_1*en_0 - s2_0*en_1)/(s2_1 - s2_0)
0137                 en_cross_B = en_0 + (BetaInverse - ri_0)*(en_1 - en_0)/(ri_1 - ri_0)
0138                 en_cross = en_cross_A if s2cross else en_cross_B
0139             else:
0140                 en_cross = np.nan
0141             pass
0142
0143             if s2_0 <= 0. and s2_1 <= 0.:
0144                 pass
0145             elif s2_0 < 0. and s2_1 > 0.:
0146                 s2integral +=  (en_1 - en_cross)*s2_1*0.5
0147             elif s2_0 >= 0. and s2_1 >= 0.:
0148                 s2integral += (en_1 - en_0)*(s2_0 + s2_1)*0.5
0149             elif s2_0 > 0. and s2_1 < 0.:
0150                 s2integral +=  (en_cross - en_0)*s2_0*0.5
0151             else:
0152                 print( " en_0 %10.5f ri_0 %10.5f s2_0 %10.5f  en_1 %10.5f ri_1 %10.5f s2_1 %10.5f " % (en_0, ri_0, s2_0, en_1, ri_1, s2_1 ))
0153                 assert 0
0154             pass
0155         pass
0156         Rfact = 369.81 / 10. #        Geant4 mm=1 cm=10
0157         NumPhotons = Rfact * charge * charge * s2integral
0158         return NumPhotons
0159
0160     def GetAverageNumberOfPhotons_s2messy(self, BetaInverse, charge=1, dump=False ):
0161         """
0162         NB see GetAverageNumberOfPhotons_s2 it gives exactly the same results as this and is simpler
0163
0164         Alternate approach doing the numerical integration directly of s2 rather than
0165         doing it on n^-2 and combining it later with the integral of the 1 and the
0166         BetaInverse*BetaInverse
0167
0168         Doing the integral of s2 avoids inconsistencies in the numerical approximations
0169         which prevents the average number of photons going negative in the region when the
0170         BetaInverse "sea level" rises to almost engulf the last rindex peak.::
0171
0172                                                   BetaInverse*BetaInverse
0173               Integral ( s2 )  = Integral(  1 - ---------------------------  )  = Integral ( 1 - c2 )
0174                                                           ri*ri
0175         """
0176         self.BetaInverse = BetaInverse
0177         ckn = self
0178         ri = ckn.ri
0179         en = ri[:,0]
0180
0181         s2 = np.zeros( (len(ri), 2), dtype=np.float64 )
0182         ct = BetaInverse/ri[:,1]
0183         s2[:,0] = ri[:,0]
0184         s2[:,1] = (1. - ct)*(1. + ct )
0185
0186         cross = ckn.FindCrossings( s2, 0. )
0187         s2integral = 0.
0188         for i in range(len(cross)//2):
0189             en0 = cross[2*i+0]
0190             en1 = cross[2*i+1]
0191             # select bins within the range
0192             s2_sel = s2[np.logical_and(s2[:,0]>=en0, s2[:,0] <= en1)]
0193
0194             # fabricate partial bins before and after the full ones
0195             # that correspond to s2 zeros
0196             fs2 = np.zeros( (2+len(s2_sel),2), dtype=np.float64 )
0197             fs2[0] = [en0, 0.]
0198             fs2[1:-1] = s2_sel
0199             fs2[-1] = [en1, 0.]
0200
0201             s2integral += np.trapz( fs2[:,1], fs2[:,0] )   # trapezoidal integration
0202         pass
0203         Rfact =  369.81   #  (eV * cm)^-1
0204         Rfact *= 0.1      # cm to mm ?  Geant4: mm = 1. cm = 10.
0205
0206         NumPhotons = Rfact * charge * charge * s2integral
0207         self.NumPhotons = NumPhotons
0208         if dump:
0209             print(" s2integral %10.4f " % (s2integral))
0210         pass
0211         return NumPhotons
0212
0213     def GetAverageNumberOfPhotons_asis(self, BetaInverse, charge=1, dump=False ):
0214         """
0215         This duplicates the results from G4Cerenkov_modified::GetAverageNumberOfPhotons
0216         including negative numbers of photons for BetaInverse close to the rindex peak.
0217
0218         Frank-Tamm formula gives number of Cerenkov photons per mm as an energy::
0219
0220                                                         BetaInverse^2
0221               N_photon  =     370.     Integral ( 1 - -----------------  )  dE
0222                                                           ri(E)^2
0223
0224         Where the integration is over regions where :    ri(E) > BetaInverse
0225         which corresponds to a real cone angle and the above bracket being positive::
0226
0227                         BetaInverse
0228               cos th = --------------   < 1
0229                            ri(E)
0230
0231         The bracket above is in fact :    1 - cos^2 th = sin^2 th  which must be +ve
0232         so getting -ve numbers of photons is clearly a bug from the numerical approximations
0233         being made.  Presumably the problem is due to the splitting of the integral into
0234         CerenkovAngleIntegral "cai" which is the cumulative integral of   1./ri(E)^2
0235         followed by linear interpolation of this in order to get the integral between
0236         crossings.
0237
0238         G4Cerenkov::
0239
0240              Rfact = 369.81/(eV * cm);
0241
0242         https://www.nevis.columbia.edu/~haas/frank_epe_course/cherenkov.ps
0243         has 370(eV.cm)^-1
0244
0245         ~/opticks_refs/nevis_cherenkov.ps
0246
0247         hc = 1240 eV nm = 1240 eV cm * 1e-7    ( nm:1e-9 cm 1e-2)
0248
0249         In [8]: 2*np.pi*1e7/(137*1240)     # fine-structure-constant 1/137 and hc = 1240 eV nm
0250         Out[8]: 369.860213514221
0251
0252         alpha/hc = 370 (eV.cm)^-1
0253
0254         See ~/opticks/examples/UseGeant4/UseGeant4.cc UseGeant4::physical_constants::
0255
0256             UseGeant4::physical_constants
0257                                                    eV 1e-06
0258                                                    cm 10
0259                                  fine_structure_const 0.00729735
0260                         one_over_fine_structure_const 137.036
0261               fine_structure_const_over_hbarc*(eV*cm) 369.81021
0262                       fine_structure_const_over_hbarc 36981020.84589
0263                             Rfact =  369.81/(eV * cm) 36981000.00000[as used by G4Cerenkov::GetAverageNumberOfPhotons]
0264                                   2*pi*1e7/(1240*137) 369.86021
0265                                                 eplus 1.00000
0266                                      electron_mass_c2 0.51099891
0267                                        proton_mass_c2 938.27201300
0268                                       neutron_mass_c2 939.56536000
0269
0270
0271         Crossing points from similar triangles::
0272
0273
0274              x - prevPM                            currentPM - prevPM
0275              ------------------------------   =     ------------------------
0276              BetaInverse - prevRI                   currentRI - prevRI
0277
0278
0279              x - prevPM =  (BetaInverse-prevRI)/(currentRI-prevRI)*(currentPM-prevPM)
0280
0281
0282
0283                              (currentPM, currentRI)
0284                                 +
0285                                /
0286                               /
0287                              /
0288                             /
0289                            /
0290                           /
0291                          *
0292                         /  (x,BetaInverse)
0293                        /
0294                       /
0295                      /
0296                     /
0297                    +
0298             (prevPM, prevRI)
0299
0300
0301         """
0302         self.BetaInverse = BetaInverse
0303         ri = self.ri
0304         cai = self.cai
0305         en = ri[:,0]
0306
0307         cross = self.FindCrossings( ri, BetaInverse )
0308         if dump:
0309             print(" cross %s " % repr(cross) )
0310         pass
0311         self.cross = cross
0312
0313         dp1 = 0.
0314         ge1 = 0.
0315
0316         for i in range(len(cross)//2):
0317             en0 = cross[2*i+0]
0318             en1 = cross[2*i+1]
0319             dp1 += en1 - en0
0320
0321             # interpolating the cai is an approximation that is the probable cause of NumPhotons
0322             # going negative for BetaInverse close to the "peak" of rindex
0323
0324             cai0 = np.interp( en0, en, cai )
0325             cai1 = np.interp( en1, en, cai )
0326             ge1 += cai1 - cai0
0327         pass
0328         Rfact =  369.81   #  (eV * cm)^-1
0329         Rfact *= 0.1      # cm to mm ?
0330
0331         NumPhotons = Rfact * charge * charge * (dp1 - ge1 * BetaInverse*BetaInverse)
0332
0333         self.dp1 = dp1
0334         self.ge1 = ge1
0335         self.NumPhotons = NumPhotons
0336
0337         if dump:
0338             print(" dp1 %10.4f ge1 %10.4f " % (dp1, ge1 ))
0339         pass
0340
0341         return NumPhotons
0342
0343     @classmethod
0344     def FindCrossings(cls, pq, pv):
0345         """
0346         :param pq: property array of shape (n,2)
0347         :param pv: scalar value
0348         :return cross: array of values where pv crosses the linear interpolated pq
0349         """
0350         assert len(pq.shape) == 2 and pq.shape[1] == 2
0351
0352         mx = pq[:,1].max()
0353         mi = pq[:,1].min()
0354
0355         cross = []
0356         if pv <= mi:
0357             cross.append( pq[0,0] )
0358             cross.append( pq[-1,0] )
0359         elif pv >= mx:
0360             pass
0361         else:
0362             if pq[0,1] >= pv:
0363                 cross.append(pq[0,0])
0364             pass
0365             assert len(pq) > 2
0366             for ii in range(1,len(pq)-1):
0367                 prevPM,prevRI = pq[ii-1]
0368                 currPM,currRI = pq[ii]
0369
0370                 down = prevRI >= pv and currRI < pv
0371                 up = prevRI < pv and currRI >= pv
0372                 if down or up:
0373                     cross.append((pv-prevRI)/(currRI-prevRI)*(currPM-prevPM) + prevPM)
0374                 pass
0375             pass
0376             if pq[-1,1] >= pv:
0377                 cross.append(pq[-1,1])
0378             pass
0379         pass
0380         assert len(cross) % 2 == 0, cross
0381         return cross
0382
0383     def test_GetAverageNumberOfPhotons(self, BetaInverse, dump=False):
0384         NumPhotons_asis = self.GetAverageNumberOfPhotons_asis(BetaInverse)
0385         NumPhotons_s2 = self.GetAverageNumberOfPhotons_s2(BetaInverse)
0386         NumPhotons_s2messy = self.GetAverageNumberOfPhotons_s2messy(BetaInverse)
0387         res = np.array([BetaInverse, NumPhotons_asis, NumPhotons_s2, NumPhotons_s2messy ])
0388         if dump:
0389             fmt = "BetaInverse %6.4f _asis %6.4f  _s2 %6.4f _s2messy %6.4f    "
0390             print( fmt % tuple(res) )
0391         pass
0392         return res
0393
0394     def scan_GetAverageNumberOfPhotons(self, x0=1., x1=2., nx=1001 ):
0395         """
0396         Creates ckn.scan comparing GetAverageNumberOfPhotons from three algorithms
0397         across the domain of BetaInverse.
0398
0399         * GetAverageNumberOfPhotons_asis
0400         * GetAverageNumberOfPhotons_s2
0401         * GetAverageNumberOfPhotons_s2messy
0402
0403         ::
0404
0405             In [12]: ckn.scan
0406             Out[12]:
0407             array([[  1.    , 293.2454, 293.2454, 293.2454],
0408                    [  1.001 , 292.7999, 292.7999, 292.7999],
0409                    [  1.002 , 292.354 , 292.354 , 292.354 ],
0410                    ...,
0411                    [  1.998 ,   0.    ,   0.    ,   0.    ],
0412                    [  1.999 ,   0.    ,   0.    ,   0.    ],
0413                    [  2.    ,   0.    ,   0.    ,   0.    ]])
0414
0415             In [13]: ckn.scan.shape
0416             Out[13]: (1001, 4)
0417
0418         """
0419         scan = np.zeros( (nx, 4), dtype=np.float64 )
0420         for i, BetaInverse in enumerate(np.linspace(x0, x1, nx )):
0421             NumPhotons_asis = self.GetAverageNumberOfPhotons_asis(BetaInverse)
0422             NumPhotons_s2 = self.GetAverageNumberOfPhotons_s2(BetaInverse, s2cross=False)
0423             NumPhotons_s2cross = self.GetAverageNumberOfPhotons_s2(BetaInverse, s2cross=True)
0424             #NumPhotons_s2messy = self.GetAverageNumberOfPhotons_s2messy(BetaInverse)
0425             scan[i] = [BetaInverse, NumPhotons_asis, NumPhotons_s2, NumPhotons_s2cross ]
0426             fmt = " bi %7.3f _asis %7.3f _s2 %7.3f _s2cross %7.3f "
0427             if i % 100 == 0: print( "%5d : %s " % (i, fmt % tuple(scan[i])) )
0428         pass
0429         self.scan = scan
0430         self.numPhotonASIS_ = lambda bi:np.interp( bi, self.scan[:,0], self.scan[:,1] )
0431         self.numPhotonS2_ = lambda bi:np.interp( bi,  self.scan[:,0], self.scan[:,2] )
0432         self.nMin = self.ri[:,1].min()
0433         self.nMax = self.ri[:,1].max()
0434
0435         d23 = scan[:,2] - scan[:,3]
0436         log.info(" d23.max %10.4f d23.min %10.4f " % (d23.max(), d23.min()))
0437
0438
0439
0440     def scan_GetAverageNumberOfPhotons_plot(self, bir=None):
0441
0442         ckn = self
0443         en = ckn.scan[:,0]
0444
0445         bi = [self.nMin, self.nMax] if bir is None else bir
0446         numPhotonMax = self.numPhotonS2_( np.linspace(bi[0], bi[1], 101) ).max()  # max in the BetaInverse range
0447
0448         extra = "%6.4f_%6.4f" % (bi[0], bi[1])
0449
0450
0451         titls = [
0452         "ana/ckn.py : scan_GetAverageNumberOfPhotons_plot %s " % extra ,
0453         "_asis goes slightly negative near rindex peak due to linear interpolation approx on top of trapezoidal integration",
0454         "_s2 avoids that by doing the integration in one pass directly on s2 (sin^2 theta) and using s2 zero crossings to improve accuracy"
0455         ]
0456
0457         title = "\n".join(titls)
0458
0459
0460         fig, ax = plt.subplots(figsize=ok.figsize)
0461         fig.suptitle(title)
0462
0463         ax.set_xlim( *bi )
0464         ax.set_ylim(  -1., numPhotonMax )
0465
0466         ax.scatter( en, ckn.scan[:,1], label="GetAverageNumberOfPhotons_asis", s=3 )
0467         ax.plot(    en, ckn.scan[:,1], label="GetAverageNumberOfPhotons_asis" )
0468
0469         ax.plot(    en, ckn.scan[:,2], label="GetAverageNumberOfPhotons_s2" )
0470         ax.scatter( en, ckn.scan[:,2], label="GetAverageNumberOfPhotons_s2", s=3 )
0471
0472         xlim = ax.get_xlim()
0473         ylim = ax.get_ylim()
0474         ax.plot( xlim, [0,0], linestyle="dotted", label="zero" )
0475
0476         for n in ckn.ri[:,1]:
0477             ax.plot( [n, n], ylim  )
0478             #ax.plot( [n, n], ylim, label="%s" % n  )
0479         pass
0480
0481         ax.legend()
0482         fig.show()
0483
0484         self.save_fig( fig, "scan_GetAverageNumberOfPhotons_plot_%s.png" % extra )
0485
0486     def save_fig(self, fig, name):
0487         path = os.path.join(self.FOLD, name)
0488         if not os.path.exists(os.path.dirname(path)):
0489             os.makedirs(os.path.dirname(path))
0490         pass
0491         log.info("save to %s " % path)
0492         fig.savefig(path)
0493         assert os.path.exists(path)
0494         self.paths.append(path)
0495
0496
0497     def scan_GetAverageNumberOfPhotons_difference_plot(self):
0498
0499         ckn = self
0500         bi = ckn.scan[:,0]
0501
0502         titls = [
0503         "ana/ckn.py : scan_GetAverageNumberOfPhotons_difference_plot : GetAverageNumberOfPhotons_s2 - GetAverageNumberOfPhotons_asis vs BetaInverse ",
0504         "refractive index values show by vertical lines",
0505         "parabolic nature of the difference because  _asis as using a linear approximation for a parabolic integral "
0506         ]
0507
0508         title = "\n".join(titls)
0509
0510         bi_range = [self.nMin, self.nMax]
0511
0512         fig, ax = plt.subplots(figsize=ok.figsize)
0513         fig.suptitle(title)
0514
0515         ax.set_xlim( *bi_range )
0516         #ax.set_ylim(  -1., numPhotonMax )
0517
0518         ax.scatter( bi, ckn.scan[:,2] - ckn.scan[:,1], s=3 )
0519         ax.plot(    bi, ckn.scan[:,2] - ckn.scan[:,1] )
0520
0521         ylim = ax.get_ylim()
0522
0523         for n in ckn.ri[:,1]:
0524             #ax.plot( [n, n], ylim, label="%s" % n  )
0525             ax.plot( [n, n], ylim )
0526         pass
0527         ax.legend()
0528         fig.show()
0529
0530         self.save_fig( fig, "scan_GetAverageNumberOfPhotons_difference_plot.png" )
0531
0532
0533     def compareOtherScans(self):
0534         """
0535         This compares the python and C++ implementations
0536         of the same double precision algorithms,
0537         agreement should be (and is) very good eg 1e-13 level
0538         """
0539         self.compare_with_QCerenkovTest()
0540         self.compare_with_G4Cerenkov_modified()
0541
0542     def compare_with_QCerenkovTest(self):
0543         """
0544         Create/update QCerenkovTest scan with::
0545
0546             qu
0547             cd tests
0548             ./QCerenkovTest.sh
0549
0550         """
0551         qck_path=os.path.expandvars("/tmp/$USER/opticks/QCerenkovTest/test_GetAverageNumberOfPhotons_s2.npy")
0552         if os.path.exists(qck_path):
0553             self.qck_scan = np.load(qck_path)
0554         else:
0555             log.error("missing %s " % qck_path)
0556         pass
0557
0558         ckn = self
0559         cf_qck_dom = np.abs( ckn.scan[:,0] - ckn.qck_scan[:,0] )
0560         cf_qck_num = np.abs( ckn.scan[:,2] - ckn.qck_scan[:,1] )
0561
0562         log.info("cf_qck_dom.min*1e6 %10.4f cf_qck_dom.max*1e6 %10.4f " % (cf_qck_dom.min()*1e6, cf_qck_dom.max()*1e6 ))
0563         log.info("cf_qck_num.min*1e6 %10.4f cf_qck_num.max*1e6 %10.4f " % (cf_qck_num.min()*1e6, cf_qck_num.max()*1e6 ))
0564
0565         assert cf_qck_dom.max() < 1e-10
0566         assert cf_qck_num.max() < 1e-10
0567
0568
0569     def compare_with_G4Cerenkov_modified(self):
0570         """
0571         Create/update scan arrays with::
0572
0573             cks
0574             ./G4Cerenkov_modifiedTest.sh
0575
0576         """
0577         cks_path="/tmp/G4Cerenkov_modifiedTest/scan_GetAverageNumberOfPhotons.npy"
0578         if os.path.exists(cks_path):
0579             self.cks_scan = np.load(cks_path)
0580         else:
0581             log.error("missing %s " % cks_path)
0582         pass
0583
0584         ckn = self
0585         cf_cks_dom = np.abs(      ckn.scan[:,0] - ckn.cks_scan[:,0] )
0586         cf_cks_num_asis = np.abs( ckn.scan[:,1] - ckn.cks_scan[:,1] )
0587         cf_cks_num_s2 = np.abs(   ckn.scan[:,2] - ckn.cks_scan[:,2] )
0588
0589         log.info("cf_cks_dom.min*1e6 %10.4f cf_cks_dom.max*1e6 %10.4f " % (cf_cks_dom.min()*1e6, cf_cks_dom.max()*1e6 ))
0590         log.info("cf_cks_num_asis.min*1e6 %10.4f cf_cks_num_asis.max*1e6 %10.4f " % (cf_cks_num_asis.min()*1e6, cf_cks_num_asis.max()*1e6 ))
0591         log.info("cf_cks_num_s2.min*1e6 %10.4f cf_cks_num_s2.max*1e6 %10.4f " % (cf_cks_num_s2.min()*1e6, cf_cks_num_s2.max()*1e6 ))
0592
0593         assert cf_cks_dom.max() < 1e-10
0594         assert cf_cks_num_asis.max() < 1e-10
0595         assert cf_cks_num_s2.max() < 1e-10
0596
0597
0598     def test_GetAverageNumberOfPhotons_plot(self, BetaInverse = 1.7):
0599         """
0600         runs test_GetAverageNumberOfPhotons for a particular BetaInverse
0601         in order to get the internals : cross, cai
0602         """
0603         ckn = self
0604         res = ckn.test_GetAverageNumberOfPhotons(BetaInverse)
0605
0606         titls = ["ana/ckn.py : test_GetAverageNumberOfPhotons_plot : %s " % str(res),
0607                  "attempt to understand how _asis manages to go negative  "
0608                  ]
0609
0610         title = "\n".join(titls)
0611
0612         cross = ckn.cross
0613         ri = ckn.ri
0614         cai = ckn.cai
0615
0616         fig, axs = plt.subplots(1, 2, figsize=ok.figsize)
0617         fig.suptitle(title)
0618
0619         ax = axs[0]
0620         ax.plot(ri[:,0], ri[:,1] , label="linear interpolation" )
0621         ax.scatter(ri[:,0], ri[:,1], label="ri" )
0622         xlim = ax.get_xlim()
0623         ylim = ax.get_ylim()
0624         ax.plot( xlim, [BetaInverse, BetaInverse], label="BetaInverse:%6.4f" % BetaInverse)
0625         for e in cross:
0626             ax.plot( [e, e], ylim, label="cross" )
0627         pass
0628
0629         ax.legend()
0630
0631         ax = axs[1]
0632         ax.plot( ri[:,0], cai, label="cai" )
0633         ylim = ax.get_ylim()
0634         for e in cross:
0635             ax.plot( [e, e], ylim, label="cross" )
0636         pass
0637         ax.legend()
0638
0639         fig.show()
0640         self.save_fig( fig, "test_GetAverageNumberOfPhotons_plot.png" )
0641
0642
0643     def load_QCerenkov_s2slv(self):
0644         """
0645         s2slv: sliver integrals across domains of BetaInverse and energy
0646
0647             In [4]: ckn.s2slv.shape
0648             Out[4]: (1001, 280)
0649
0650         cs2slv: is cumulative sum of the above sliver integrals with the same shape
0651
0652             In [5]: ckn.cs2slv.shape
0653             Out[5]: (1001, 280)
0654
0655         s2slv_sum: sum over energy domain of s2slv
0656
0657             In [8]: ckn.s2slv_sum.shape
0658             Out[8]: (1001,)
0659
0660         cs2slv_last: last of the (energy axis) cumulative sum values : this very closely matches s2slv_sum
0661
0662             In [9]: ckn.cs2slv_last.shape
0663             Out[9]: (1001,)
0664
0665         c2slv_norm: energy axis cumsum divided by the last : making this the CDF
0666
0667             In [10]: ckn.cs2slv_norm.shape
0668             Out[10]: (1001, 280)
0669
0670             Notice lots of nan, for CK disallowed BetaInverse
0671
0672         """
0673         s2slv_path = "/tmp/blyth/opticks/QCerenkovTest/test_getS2SliverIntegrals_many.npy"
0674         log.info("load %s " % s2slv_path )
0675         s2slv = np.load(s2slv_path) if os.path.exists(s2slv_path) else None
0676         globals()["ss"] = s2slv
0677
0678         cs2slv_path = "/tmp/blyth/opticks/QCerenkovTest/test_getS2SliverIntegrals_many_cs2slv.npy"
0679         log.info("load %s " % cs2slv_path )
0680         cs2slv = np.load(cs2slv_path) if os.path.exists(cs2slv_path) else None
0681         cs2slv_last = cs2slv[:,-1]
0682
0683         cs2slv_norm_path = "/tmp/blyth/opticks/QCerenkovTest/test_getS2SliverIntegrals_many_cs2slv_norm.npy"
0684         log.info("load %s " % cs2slv_norm_path )
0685         cs2slv_norm = np.load(cs2slv_norm_path) if os.path.exists(cs2slv_norm_path) else None
0686
0687         globals()["cs"] = cs2slv
0688         globals()["csl"] = cs2slv_last
0689         globals()["csn"] = cs2slv_norm
0690
0691         self.s2slv = s2slv
0692         self.cs2slv = cs2slv
0693         self.cs2slv_last = cs2slv_last
0694         self.cs2slv_norm = cs2slv_norm
0695
0696         s2slv_sum = s2slv.sum(axis=1)
0697         self.s2slv_sum = s2slv_sum
0698
0699         sum_vs_last = np.abs( s2slv_sum - cs2slv_last )
0700         sum_vs_last_mx = sum_vs_last.max()
0701         log.info( "sum_vs_last_mx %s " % sum_vs_last_mx )
0702         assert sum_vs_last_mx < 1e-10
0703
0704
0705
0706     def load_QCerenkov_s2slv_plot(self):
0707         """
0708
0709         """
0710         titls = [
0711         "ana/ckn.py : load_QCerenkov_s2slv_plot : compare sum of s2slv from many sliver bins to the rindex big bin result",
0712         "deviation of up to 0.4 of a photon between the sum of sliver integrals and the big bin integrals is as yet unexplained"
0713         ]
0714         title = "\n".join(titls)
0715
0716         fig, ax = plt.subplots(figsize=ok.figsize)
0717         fig.suptitle(title)
0718
0719         ax.plot( self.qck_scan[:,0], self.qck_scan[:,1], label="qck_scan" )
0720         ax.plot( self.qck_scan[:,0], self.s2slv_sum,     label="s2slv_sum" )
0721         ax.plot( self.qck_scan[:,0], self.cs2slv_last,   label="cs2slv_last" )
0722         ax.legend()
0723
0724         axr = ax.twinx()
0725         axr.plot( self.qck_scan[:,0], self.s2slv_sum - self.qck_scan[:,1], label="diff", linestyle="dotted" )
0726         axr.legend(loc='lower right')
0727
0728         fig.show()
0729
0730         self.ax = ax
0731
0732
0733 if __name__ == '__main__':
0734
0735     ok = opticks_main()
0736     ckn = CKN()
0737
0738     ckn.scan_GetAverageNumberOfPhotons()
0739
0740     ckn.scan_GetAverageNumberOfPhotons_plot(bir=[1,2])
0741     ckn.scan_GetAverageNumberOfPhotons_plot()
0742     ckn.scan_GetAverageNumberOfPhotons_plot(bir=[1.7,1.8])
0743     ckn.scan_GetAverageNumberOfPhotons_difference_plot()
0744
0745     ckn.compareOtherScans()
0746
0747     ckn.test_GetAverageNumberOfPhotons_plot(BetaInverse=1.788)
0748
0749     ckn.load_QCerenkov_s2slv()
0750     ckn.load_QCerenkov_s2slv_plot()
0751
0752     print("\n".join(ckn.paths))
0753
0754
0755
0756