eic-opticks/ana/ckphotons.py

0001 #!/usr/bin/env python
0002
0003 import math
0004 import numpy as np
0005
0006 from scipy.interpolate import interp1d
0007 from opticks.ana.mlib import GMaterialLib
0008 mlib = GMaterialLib()
0009
0010 import matplotlib.pyplot as plt
0011
0012
0013 class G4Cerenkov(object):
0014     """
0015     391 void G4Cerenkov::BuildThePhysicsTable()
0016     392 {
0017     ...
0018     424             // Retrieve the first refraction index in vector
0019     425             // of (photon energy, refraction index) pairs
0020     426
0021     427             G4double currentRI = (*theRefractionIndexVector)[0];
0022     428
0023     429             if (currentRI > 1.0) {
0024     430
0025     431                // Create first (photon energy, Cerenkov Integral)
0026     432                // pair
0027     433
0028     434                G4double currentPM = theRefractionIndexVector->Energy(0);
0029     435                G4double currentCAI = 0.0;
0030     436
0031     437                aPhysicsOrderedFreeVector->InsertValues(currentPM , currentCAI);
0032     438
0033     439                // Set previous values to current ones prior to loop
0034     440
0035     441                G4double prevPM  = currentPM;
0036     442                G4double prevCAI = currentCAI;
0037     443                G4double prevRI  = currentRI;
0038     444
0039     445                // loop over all (photon energy, refraction index)
0040     446                // pairs stored for this material
0041     447
0042     448                for (size_t ii = 1;
0043     449                            ii < theRefractionIndexVector->GetVectorLength();
0044     450                            ++ii) {
0045     451                    currentRI = (*theRefractionIndexVector)[ii];
0046     452                    currentPM = theRefractionIndexVector->Energy(ii);
0047     453
0048     454                    currentCAI = 0.5*(1.0/(prevRI*prevRI) +
0049     455                                      1.0/(currentRI*currentRI));
0050     456
0051     457                    currentCAI = prevCAI + (currentPM - prevPM) * currentCAI;
0052     458
0053     459                    aPhysicsOrderedFreeVector->
0054     460                                          InsertValues(currentPM, currentCAI);
0055     461
0056     462                    prevPM  = currentPM;
0057     463                    prevCAI = currentCAI;
0058     464                    prevRI  = currentRI;
0059     465                }
0060     466
0061     467             }
0062     ...
0063     """
0064     @classmethod
0065     def CerenkovAngleIntegral(cls, ev, ri):
0066         """
0067         :param ev: np.ndarray domain of energy (eV) in ascending order
0068         :param ri: np.ndarry values of RINDEX on the domain
0069         :return cai:  np.ndarray numerical 1/ri^2 energy integral starting at zero
0070         """
0071         assert len(ev) == len(ri)
0072         x = ev
0073         y = 1./(ri*ri)
0074         ymid = (y[:-1]+y[1:])/2.                     # looses entry as needs pair
0075         xdif = np.diff(x)                            # also looses entry
0076         cai = np.zeros( len(ri), dtype=np.float32 )  # gains entry from first zero
0077         np.cumsum(ymid*xdif, out=cai[1:] )
0078         return cai
0079     pass
0080
0081
0082
0083 if __name__ == '__main__':
0084
0085     ri0 = mlib("Water.RINDEX").copy()
0086
0087     # energy(eV) and refractive index in ascending energy
0088     ev = mlib.ev[::-1]
0089     ri = ri0[::-1]
0090
0091     fig, ax = plt.subplots(1)
0092     ax.plot( ev, ri, "o" )
0093     fig.show()
0094
0095
0096 if 0:
0097
0098     cai = G4Cerenkov.CerenkovAngleIntegral(ev, ri)
0099
0100     eV = 1.0e-06
0101     cm = 10.
0102     Rfact = 369.81/(eV * cm)
0103
0104     MeV = 1.
0105     GeV = 1000.
0106
0107     mass_c2 = {}
0108     mass_c2["electron"] = 0.510998910 * MeV
0109     mass_c2["proton"] = 938.272013 * MeV
0110     mass_c2["muon"] = 105.659 * MeV
0111
0112     eplus = 1.0
0113     charge = eplus
0114
0115     e_total = 200.*GeV
0116     e_rest = mass_c2["muon"]
0117     gamma = e_total/e_rest
0118     beta = math.sqrt( 1. - 1./(gamma*gamma) )
0119     BetaInverse = 1./beta
0120
0121     CAImax = cai.max()
0122
0123     Pmin = ev[0]
0124     Pmax = ev[-1]
0125
0126     nMin = ri.min()
0127     nMax = ri.max()
0128
0129     print("%-20s : %s " % ("BetaInverse", BetaInverse))
0130     print("%-20s : %s " % ("nMin", nMin))
0131     print("%-20s : %s " % ("nMax", nMax))
0132
0133     if nMax < BetaInverse:
0134         dp = 0.0
0135         ge = 0.0
0136     elif nMin > BetaInverse:
0137         dp = Pmax - Pmin
0138         ge = CAImax
0139     else:
0140         """
0141         //  If n(Pmin) < 1/Beta, and n(Pmax) >= 1/Beta, then
0142         // we need to find a P such that the value of n(P) == 1/Beta.
0143         // Interpolation is performed by the GetEnergy() and
0144         // Value() methods of the G4MaterialPropertiesTable and
0145         // the GetValue() method of G4PhysicsVector.
0146         //Pmin = Rindex->GetEnergy(BetaInverse);
0147         """
0148         assert np.all( np.diff(ri) > 0. )   # usually ri is not monotonic
0149         Pmin = np.interp( BetaInverse, ri, ev )
0150         dp = Pmax - Pmin
0151
0152         """
0153         // need boolean for current implementation of G4PhysicsVector
0154         // ==> being phased out
0155         G4bool isOutRange;
0156         G4double CAImin = CerenkovAngleIntegrals->GetValue(Pmin, isOutRange);
0157         ge = CAImax - CAImin;
0158         """
0159         assert np.all( np.diff(ev) > 0. )
0160         CAImin = np.interp( Pmin, ev, cai )
0161         ge = CAImax - CAImin
0162     pass
0163     NumPhotons = Rfact * charge/eplus * charge/eplus * (dp - ge * BetaInverse*BetaInverse)
0164
0165
0166
0167
0168
0169