Commit 2668bc4c authored by Yakubov, Sergey's avatar Yakubov, Sergey
Browse files

add lib folder

parent 57ae1048

.gitignore

+0 −2
Original line number Diff line number Diff line
@@ -12,8 +12,6 @@ build/ 
develop-eggs/

dist/

eggs/

lib/

lib64/

parts/

sdist/

var/

lib/ISAW_PythonSources/Lib/SNSlibs.py

0 → 100644
+376 −0
Original line number Diff line number Diff line 


from DataSetTools.operator import Operator

from DataSetTools.operator.Generic.Load import GenericLoad

from gov.anl.ipns.Parameters import * # parameters

from gov.anl.ipns.Util.SpecialStrings import ErrorString

from Operators.TOF_Diffractometer import *

from DataSetTools.operator.DataSet.EditList import *

from Operators.Special import *

from Operators.TOF_SCD import *

from Command import ScriptUtil

from NexIO.NexApi import *;

from NexIO.Util import ConvertDataTypes;

from java.lang import *



def getIsawHome():

        

        Res = System.getProperty("ISAW_HOME")

       

        if not Res.endswith("/"):

           Res +=  '/'

        return Res



def getDefaultDetCalFile( instr):

         Res = getIsawHome()+"InstrumentInfo/SNS/"+instr+"/"+instr+".DetCal"

         return Res



def getDefaultBankFile( instr):

         Res = getIsawHome()+"InstrumentInfo/SNS/"+instr+"/"+instr+"_bank.xml"

         return Res



def getDefaultMapFile( instr):

         Res = getIsawHome()+"InstrumentInfo/SNS/"+instr+"/"+instr+"_TS.dat"

         return Res



def getDefaultVanPeakFile( ):

         Res = getIsawHome()+"Databases/VanadiumPeaks.dat"

         return Res







'''def getCalibPG(instr, type):

        cmd = "findcalib -i"+instr+" --listall " + type

        print cmd

        from java.lang import Runtime

        proc = Runtime.getRuntime().exec(cmd)

        print proc

        import ExtTools.monq.stuff.Exec

        thing = ExtTools.monq.stuff.Exec(proc)

        if not thing.done():

            pass

        print thing.getOutputText()

        return None # should make a system call to findcalib

'''

def sumGroups( inds, outds, groups, newId, position):

        data = inds.getData_entry_with_id(groups[0])

        if data is None:

            raise "Something went wrong with group id = %d" % groups[0]

        for i in groups[1:]:

            data = data.add(inds.getData_entry_with_id(i))

            if data is None:

                raise "Something went wrong with group id = %d" % i

        data.setGroup_ID(newId)

        data.setAttribute(position)

        outds.addData_entry(data)



def getPosAttr( bragg, L2):

        from gov.anl.ipns.MathTools.Geometry import Vector3D, DetectorPosition

        import math

        pos = [L2*math.cos(bragg*math.pi/180.),

               L2*math.sin(bragg*math.pi/180.),

               0., 0.] # last one is the weight

        pos = Vector3D(pos)

        pos = DetectorPosition(pos)



        from DataSetTools.dataset import DetPosAttribute

        attr = DetPosAttribute("Effective Position", pos)

        

        return attr



def getRunDir(instr, runnumber):

        cmd = "findnexus -i"+instr+" -A %d" % runnumber

        import os

        f = os.popen(cmd)

        nxs = f.readline()

        nxs = nxs.strip()

        if not os.path.exists(nxs):

            raise Error("Failed to find run \"%d\": %s"  % (runnumber, nxs))

        temp = os.path.split(nxs)[0] # NeXus directory

        temp = os.path.split(temp)[0] # run directory

        return (nxs, os.path.join(temp, "preNeXus"))



def getRunStuff(instr, runnumber):

        import os

       

        # get the nexus file and prenexus directory

        (nxs, prenxs) = getRunDir(instr, runnumber)

        # print ["nxs,prenxs",nxs,prenxs]

        # determine the event file

        event = "%s_%d_neutron_event.dat" % (instr, runnumber)

        event = os.path.join(prenxs, event)

        #print ["event",event]

#        if not os.path.exists(event):

#            raise Error, "%s does not exist" % event



        # determine the cvinfo file

        cvinfo = "%s_%d_cvinfo.xml" % (instr, runnumber)

        cvinfo = os.path.join(prenxs, cvinfo)

        if not os.path.exists(cvinfo):

            raise Error("%s does not exist" % cvinfo)



        # get the protoncharge

        pcharge = 1.

        cvinfo = open(cvinfo, "r")

        temp = cvinfo.readline()

        while not "protoncharge" in temp:

            temp = cvinfo.readline()

        if len(temp) > 0:

            start = temp.index("value=\"") + len("value=\"")

            stop = temp.index("\"", start + 1)

            pcharge = float(temp[start:stop])

        cvinfo.close()



        return (event, pcharge)



def fixUnits( dataset):

        """This needs to be done if the label and units were interchanged"""

        units = dataset.getX_units()

        label = dataset.getX_label()

        if units == "Time-of-flight":

            dataset.setX_units(label)

            dataset.setX_label(units)



def removeZeros( dataset):

        ids = list(range(dataset.getNum_entries()))

        for i in ids:

            spectrum = dataset.getData_entry(i)

            y = spectrum.getY_values()

            dy = spectrum.getErrors()

            for j in range(len(y)):

                if y[j] == 0.:

                    y[j] = 1.

                    dy[j] = 1.

            #dataset.replaceData_entry(spectrum, i)



def send( ds, showPlots, sendData,IOBS):

        if ds is None:

            return

        from Command import ScriptUtil

        if showPlots:

            ScriptUtil.display(ds)

        if sendData and IOBS is not None:

            ScriptUtil.send(ds, IOBS)



def EventD_space(sendData,showData,IOBS,instr,runnum,DetCalFile,BankFile,

             MapFile,firstEvent,NumEvents,d_min,d_max,log_param,scale):

        

        (eventFile, pcharge) = getRunStuff(instr,runnum)

        

        pcharge = getProtonsCharge( instr, runnum)

        pcharge = pcharge/scale

        print("scale = %.2f" % pcharge)





        log_param = log_param * d_min # log param

        useLogBinning = 1

        nbins =0

        useD_spaceMapFile = 0

        d_spaceMapFile =""

        # do the initial load and time focus

        import os



#        if not os.path.exists(event):

#            raise Error, "%s does not exist" % event

        

        args= [eventFile,DetCalFile,BankFile,MapFile,firstEvent,

                                  NumEvents,d_min,d_max,useLogBinning,log_param,nbins,

                                  0,"",0,"",0,0]

        if  os.path.exists(eventFile):

             panel_ds = Util.Make_d_DataSet( eventFile,DetCalFile,BankFile,MapFile,firstEvent,NumEvents,d_min,d_max,useLogBinning,log_param,nbins, 0,"",0,"",0,0)

        else:

             L = eventFile.split('_')

             L[len(L)-2]="neutron0"

             eventFile1 = ""

             for i in range(0,len(L)-1):

                eventFile1 += L[i]+'_'

             eventFile1 += L[len(L)-1]

             panel_ds1 = Util.Make_d_DataSet( eventFile1,DetCalFile,BankFile,MapFile,firstEvent,NumEvents,d_min,d_max,useLogBinning,log_param,nbins, 0,"",0,"",0,0)

             eventFile2 =""

             L[len(L)-2]="neutron1"

             for i in range(0,len(L)-1):

                eventFile2 += L[i]+'_'

             eventFile2 += L[len(L)-1]

             panel_ds2 = Util.Make_d_DataSet( eventFile2,DetCalFile,BankFile,MapFile,firstEvent,NumEvents,d_min,d_max,useLogBinning,log_param,nbins, 0,"",0,"",0,0)

             panel_ds = DataSetMerge(panel_ds1,panel_ds2).getResult()     



        

       

        send(panel_ds, showData, sendData,IOBS)



        panel_ds.setSqrtErrorsAtLeast_1()



        return panel_ds



def EventD_space2GSAS(sendData,showData,IOBS,instr,runnum,DetCalFile,BankFile,

             MapFile,firstEvent,NumEvents,d_min,d_max,log_param,scale):

        

        (eventFile, pcharge) = getRunStuff(instr,runnum)

        

        pcharge = getProtonsCharge( instr, runnum)

        pcharge = pcharge/scale

        print("scale = %.2f" % pcharge)





        log_param = log_param * d_min # log param

        useLogBinning = 1

        nbins =0

        useD_spaceMapFile = 0

        d_spaceMapFile =""

        # do the initial load and time focus

        import os



#        if not os.path.exists(event):

#            raise Error, "%s does not exist" % event

        

        args= [eventFile,DetCalFile,BankFile,MapFile,firstEvent,

                                  NumEvents,d_min,d_max,useLogBinning,log_param,nbins,

                                  0,"",0,"",0,0]

        if  os.path.exists(eventFile):

             panel_ds = Util.Make_d_DataSet( eventFile,DetCalFile,BankFile,MapFile,firstEvent,NumEvents,d_min,d_max,useLogBinning,log_param,nbins, 0,"",0,"",0,0)

        else:

             L = eventFile.split('_')

             L[len(L)-2]="neutron0"

             eventFile1 = ""

             for i in range(0,len(L)-1):

                eventFile1 += L[i]+'_'

             eventFile1 += L[len(L)-1]

             panel_ds1 = Util.Make_d_DataSet( eventFile1,DetCalFile,BankFile,MapFile,firstEvent,NumEvents,d_min,d_max,useLogBinning,log_param,nbins, 0,"",0,"",0,0)

             eventFile2 =""

             L[len(L)-2]="neutron1"

             for i in range(0,len(L)-1):

                eventFile2 += L[i]+'_'

             eventFile2 += L[len(L)-1]

             panel_ds2 = Util.Make_d_DataSet( eventFile2,DetCalFile,BankFile,MapFile,firstEvent,NumEvents,d_min,d_max,useLogBinning,log_param,nbins, 0,"",0,"",0,0)

             panel_ds = DataSetMerge(panel_ds1,panel_ds2).getResult()     



        

       

        send(panel_ds, showData, sendData,IOBS)



        panel_ds.setSqrtErrorsAtLeast_1()

        

        

        # sum spectra together

        gsas_ds = panel_ds.clone()

        

        # convert back to tof

        op = gsas_ds.getOperator("ToTof")

        if op is None:

            op = gsas_ds.getOperator("Convert d-Spacing to TOF")

        gsas_tof = op.getResult()

        if isinstance(gsas_tof, ErrorString): # error occured

            return gsas_tof

        title = gsas_tof.getTitle().replace("_d-spacing", " tof")

        gsas_tof.setTitle(title)

        send(gsas_tof, showData, sendData,IOBS)



        # normalize the data

        op = gsas_tof.getOperator("Divide by Scalar")

        scalar = op.getParameter(0)

        import java.lang.Float

        scalar.setValue(java.lang.Float(pcharge))

        op.getResult()

        send(gsas_tof, showData, sendData,IOBS)

       

        fixUnits(gsas_tof)

        removeZeros(gsas_tof)



        send( gsas_tof, showData,sendData,IOBS)

        return gsas_tof



def run():

    print(EventD_space2GSAS( 1,0,IOBS,"PG3",666,"/SNS/users/pf9/NEW_PG3.DetCal","/SNS/PG3/2010_2_11_CAL/calibrations/PG3_bank_2010_04_22.xml","/SNS/PG3/2009_2_11A_CAL/calibrations/PG3_TS_2009_04_17.dat",0,1E12,.2,5,2E-4,9.99E12))



def FixVanadium( sendData,showData,IOBS,instr,runnum,DetCalFile,BankFile,

             MapFile,firstEvent,NumEvents,d_min,d_max,log_param,scale,BadPeakFile,

             PeakWidth_bad,PeakInterval_bad,NChanAv_bad,CutOffFilter,OrderFilter

               ):

        # whether or not to send all datasets to the tree

       



        

        print(firstEvent)

        ds =EventD_space2GSAS(0,showData,IOBS,instr,runnum,DetCalFile,BankFile,

             MapFile,firstEvent,NumEvents,d_min,d_max,log_param,scale)

       

      

        send(ds.clone(), showData, sendData,IOBS)

        print(["ds?", ds.getNum_entries()])

        RemovePeaks_Calc.RemovePeaks_tof(ds,BadPeakFile,PeakWidth_bad, PeakInterval_bad,NChanAv_bad)



        

        Res = LowPassFilterDS(ds, CutOffFilter, OrderFilter).getResult()



        if isinstance( Res, ErrorString):

           return Res



        send(ds, showData, sendData,IOBS)

       

        return ds

        

def run5():

     print(FixVanadium(1,0,IOBS,"PG3",539,"/SNS/users/pf9/NEW_PG3.DetCal",

       "/SNS/PG3/2010_2_11_CAL/calibrations/PG3_bank_2010_03_11.xml",

        "/SNS/PG3/2010_2_11A_CAL/calibrations/PG3_TS_2009_04_17.dat",

        0,1E12,.2,10,2E-4,9.999E10,"/SNS/users/ehx/SNS_ISAW/Databases/VanadiumPeaks.dat",.005,1.9,10,.02,2))



def rotateDetectors( instr, runnum, DetCalFile1):

     

      if DetCalFile1 is None:

          DetCalFile1 = getDefaultDetCalFile( instr)

      

      import os

      (nxs, prenxs) = getRunDir(instr, runnum)

      cvinfo = "%s_%d_cvinfo.xml" % (instr, runnum)

      cvinfo = os.path.join(prenxs, cvinfo)



      print("CVinfo=",cvinfo,runnum,instr)

      if not os.path.exists(cvinfo):

         raise Error("%s does not exist" % cvinfo)

           

      V = Util.getRotationAngles( cvinfo)

      if V is None or V.size()< 1:

          print("No information to rotate the Detectors")

          return DetCalFile1



      ang1 = V.firstElement().firstElement()

      ang2 = V.elementAt(1).firstElement()

    

      

      DetCal1 = System.getProperty("user.home")

      if not DetCal1.endswith("/"):

          DetCal1 +='/'

      Save = DetCal1

      Save= DetCal1

      DetCal1 += "ISAW/tmp/dummy.DetCal"

      DetCal2 = Save+"ISAW/tmp/dummy2.DetCal"



      General_Utils.RotateDetectors( DetCalFile1,14,DetCal1, ang1,0,0)

      try:

         General_Utils.RotateDetectors( DetCal1,5,DetCal2, ang2,0,0, ScriptUtil.ToVec([1,2,3,4,5,6,7,8,9]))

      except:

         return DetCal1   





      return DetCal2



def getProtonsCharge( instr, runnum):

      import os

      (nxs, prenxs) = getRunDir(instr, runnum)

     

      node = NexNode( nxs)     

      for i in range(0,node.getNChildNodes()):

       

         n = node.getChildNode( i)       

        

         if n.getNodeClass()=="NXentry":

            

            try:

               V = n.getChildNode("proton_charge").getNodeValue()

              

               F = ConvertDataTypes.floatValue(V)

               return Float(F).doubleValue()

            except:

               pass



      return Double.NaN

lib/ISAW_PythonSources/Lib/SNSlibs.py.bak

0 → 100644
+376 −0

File added.

Preview size limit exceeded, changes collapsed.

lib/ISAW_PythonSources/Lib/absor_V_rod.py

0 → 100644
+100 −0
Original line number Diff line number Diff line 
#--------------------------------------------------------

#               function absor_V_rod

#--------------------------------------------------------

# Function to calculate the absorption and of a vanadium rod.

#--------------------------------------------------------

#   A.J. Schultz,   September, 2010

#--------------------------------------------------------

# Subroutine to calculate the absorption correction

# for a vanadium rod. Based on values for cylinders in:

# 

# C. W. Dwiggins, Jr., Acta Cryst. A31, 146 (1975).

# 

# In this paper, A is the transmission and A* = 1/A is

# the absorption correction.

# 

# For each of the 19 theta values in Dwiggins (theta = 0.0 to 90.0

# in steps of 5.0 deg.), the ASTAR values vs.muR were fit to a fourth

# order polynomial in Excel. These values are given below in the

# data statement. (For a sphere, third order polynomials were

# sufficient, but not for the cylinder.)



from math import *



def absor_V_rod( angleH, angleV, wl):

    "Returns the absorption correction for a vanadium cylinder."

    

    # pc = fourth order polynomial coefficients

    pc = [ [ 1.000,   1.6516,   1.6251,   0.2971,   0.5155 ],

           [ 1.000,   1.7483,   1.2967,   0.6254,   0.3947 ],

           [ 1.000,   1.9176,   0.6903,   1.2585,   0.1359 ],

           [ 1.000,   1.9995,   0.3439,   1.6594,  -0.0884 ],

           [ 1.000,   1.9627,   0.3897,   1.6581,  -0.2035 ],

           [ 1.000,   1.8675,   0.6536,   1.3886,  -0.2274 ],

           [ 1.000,   1.7628,   0.9709,   1.0092,  -0.1983 ],

           [ 1.000,   1.6821,   1.2256,   0.6416,  -0.1503 ],

           [ 1.000,   1.6336,   1.3794,   0.3412,  -0.1026 ],

           [ 1.000,   1.6114,   1.4430,   0.1151,  -0.0620 ],

           [ 1.000,   1.6075,   1.4387,  -0.0451,  -0.0303 ],

           [ 1.000,   1.6146,   1.3900,  -0.1534,  -0.0067 ],

           [ 1.000,   1.6278,   1.3150,  -0.2219,   0.0097 ],

           [ 1.000,   1.6406,   1.2405,  -0.2719,   0.0230 ],

           [ 1.000,   1.6566,   1.1508,  -0.2901,   0.0293 ],

           [ 1.000,   1.7516,   0.8713,  -0.1695,   0.0086 ],

           [ 1.000,   1.6771,   1.0277,  -0.3139,   0.0380 ],

           [ 1.000,   1.6826,   0.9942,  -0.3185,   0.0399 ],

           [ 1.000,   1.6851,   0.9814,  -0.3190,   0.0403 ] ]



    # From Structure of Metals by Barrett and Massalksi:

    #                vanadium is b.c.c, a = 3.0282

    # V = 27.769, Z = 2

    # From lin_abs_coef in ISAW:

    smu = 0.367  # linear absorption coeff. for total scattering in cm_1

    amu = 0.366  # linear absorption coeff. for true absorption at 1.8 A in cm^-1

    radius = 0.407  # radius of the vanadium rod used for TOPAZ

    

    mu = smu + (amu/1.8)*wl



    muR = mu*radius

    

    theta = (angleH*180.0/pi)/2.0   # theta is the theta angle in the horizontal plane



    # ! Using the polymial coefficients, calulate ASTAR (= 1/transmission) at

    # ! theta values below and above the actual theta value.



    i = int(theta/5.0)

    astar1 = pc[i][0] + pc[i][1]*muR + pc[i][2]*muR**2 + pc[i][3]*muR**3



    i = i+1

    astar2 = pc[i][0] + pc[i][1]*muR + pc[i][2]*muR**2 + pc[i][3]*muR**3



    # !	Do a linear interpolation between theta values.



    frac = (theta%5.0)/5.0



    astar = astar1*(1-frac) + astar2*frac	# astar is the correction



    trans1 = 1.0/astar	                        # trans is the transmission

                                                # trans = exp(-mu*tbar)

	

    # !	Calculate TBAR as defined by Coppens.



    tbar1 = -log(trans1)/mu



    # Calculate total path length and transmission for scattered

    # beam out of the horizontal plane.

    

    tbar2 = tbar1 / cos( angleV )  # path length

    trans2 = exp( -mu * tbar2 )    # transmission

    

    return trans2



# test

    

# CenterX = -23.2175

# CenterY = 0.00

# CenterZ = 31.962

# wl = 2.0



# transmission = absor_V_rod( CenterX, CenterY, CenterZ, wl )

# print transmission

lib/ISAW_PythonSources/Lib/absor_sphere.py

0 → 100644
+97 −0
Original line number Diff line number Diff line 
#--------------------------------------------------------

#               function absor_sphere

#--------------------------------------------------------

# Function to calculate the absorption and  of a sherical

# crystal.

#--------------------------------------------------------

# Jython version:

#   A.J. Schultz,   November, 2009

#--------------------------------------------------------

# Comments from the Fortran source code in anvredSNS.f:

# ! Subroutine to calculate a spherical absorption correction

# ! and tbar. Based on values in:

# !

# ! C. W. Dwiggins, Jr., Acta Cryst. A31, 395 (1975).

# !

# ! In this paper, A is the transmission and A* = 1/A is

# ! the absorption correction.

#

# ! Input are the smu (scattering) and amu (absorption at 1.8 Ang.)

# ! linear absorption coefficients, the radius R of the sample

# ! the theta angle and wavelength.

# ! The absorption (absn) and tbar are returned.

#

# ! A. J. Schultz, June, 2008

#

#   real mu, muR    !mu is the linear absorption coefficient,

#                   !R is the radius of the spherical sample.

#

# ! For each of the 19 theta values in Dwiggins (theta = 0.0 to 90.0

# ! in steps of 5.0 deg.), the ASTAR values vs.muR were fit to a third

# ! order polynomial in Excel. These values are given below in the

# ! data statement.



from math import *



def absor_sphere(smu, amu, radius, twoth, wl):

    "Returns the absorption correction for a Bragg peak."

    

    

    # pc = third order polynomial coefficients

    pc = [ [ 1.0000,  1.9368,  0.0145,  1.1386 ],

           [ 1.0000,  1.8653,  0.1596,  1.0604 ],

           [ 1.0000,  1.6908,  0.5175,  0.8598 ],

           [ 1.0000,  1.4981,  0.9237,  0.6111 ],

           [ 1.0000,  1.3532,  1.2436,  0.3798 ],

           [ 1.0000,  1.2746,  1.4308,  0.1962 ],

           [ 1.0000,  1.2530,  1.4944,  0.0652 ],

           [ 1.0000,  1.2714,  1.4635, -0.0198 ],

           [ 1.0000,  1.3093,  1.3770, -0.0716 ],

           [ 1.0000,  1.3559,  1.2585, -0.0993 ],

           [ 1.0000,  1.4019,  1.1297, -0.1176 ],

           [ 1.0000,  1.4434,  1.0026, -0.1153 ],

           [ 1.0000,  1.4794,  0.8828, -0.1125 ],

           [ 1.0000,  1.5088,  0.7768, -0.1073 ],

           [ 1.0000,  1.5317,  0.6875, -0.1016 ],

           [ 1.0000,  1.5489,  0.6159, -0.0962 ],

           [ 1.0000,  1.5608,  0.5637, -0.0922 ],

           [ 1.0000,  1.5677,  0.5320, -0.0898 ],

           [ 1.0000,  1.5700,  0.5216, -0.0892 ] ]



    mu = smu + (amu/1.8)*wl



    muR = mu*radius

    

    theta = (twoth*180.0/pi)/2.0



# ! Using the polymial coefficients, calulate ASTAR (= 1/transmission) at

# ! theta values below and above the actual theta value.



    i = int(theta/5.0)

    astar1 = pc[i][0] + pc[i][1]*muR + pc[i][2]*muR**2 + pc[i][3]*muR**3



    i = i+1

    astar2 = pc[i][0] + pc[i][1]*muR + pc[i][2]*muR**2 + pc[i][3]*muR**3



# ! Do a linear interpolation between theta values.



    frac = (theta%5.0)/5.0



    astar = astar1*(1-frac) + astar2*frac	# astar is the correction



    trans = 1.0/astar                          # trans is the transmission

                                               # trans = exp(-mu*tbar)



# ! Calculate TBAR as defined by Coppens.



    if mu == 0.0:

        tbar = 0.0

    else:

        tbar = -log(trans)/mu





    return trans, tbar