Getting Calcite constants out of TDS, a worked out example

On the ftp site where this documentation is, you can find first the input file for peaks extraction extraction.par here how it looks

images_prefix             =  "data/Calcite_170K_1_0"

# Let the following images_prefix_alt to None and dont set it ( default None)
# It is used for reviewing with comparison of same position locations
# versus an alternative dataset
#  images_prefix_alt         =  "G19K45/set1_0001p_"
images_postfix        =  ".cbf"

## peaks are searched, having a radius peak_size
#  and with maximum > threshold_abs and > image_max*threshold_rel
threshold_abs=40000
threshold_rel=0.00
peak_size = 2
spot_view_size = 10

# Filtering, OPTIONAL or set it to None for no mask filtering of bad pixels
# filter_file           =  None  
filter_file           =  "data/filter.edf"  

##  medianize can be 1 or 0 ( True or False)
##  The peak is always searched into a stack
##  of spot_view_size + 1 + spot_view_size slices
##  When the option is activated each slice will be beforehand
##  subtracted with the median of the  slices going from the spot_view_size^th preceeding
##  to the spot_view_size^th following. Spot_view_size should be greaer than peak_size
medianize = 1


## Let this advanced option as it is for now, or dive into the code.
## The rationale here would be that in the first step on works
## locally on spot_view_size + 1 + spot_view_size slices,
##  in the second step one can work globally on the found peaks
globalize = (0,0)

## The file where  found peaks are written 
## or, when reviewing, read. 
file_maxima = "maxima.p"

### if yes, a graphical window allows for the selection of the mask
### Closing the window will propose to write on newmask.edf
drawmask = False

## Allows to have a simultaneous and animated view of the blobs.
## Animated thumnails.
## They can be selected of deselected. When clicking on the
## thumbnails a larger 3D explorer is charged with
## a (review_size+1+review_size )^3 stack
review = False
review_size = 50
that you use on the data that you can download(3G!!!) here
calcite.tgz.

You have to create a directory test and you write there extraction.par and you untar there calcite.tgz. This will create the image files in subdirectory data of test.

Do first a peak detection run, that will extract approximate peak positions which will be used for geometry alignement

tds2el_extraction_maxima_vX  extraction.par

where vX is the version number (you can set it in __init.py before installation).

This generates file maxima.p as requested. We are going to use these peak positions to align the sample. You can check that all the peaks are meaningful for an alignement by setting review to True and rerun, but we are not discarding any for the moment(if you want you can already view the peaks by setting review to True but dont discard them unless there is a good reason because the alignement of the sample runs better with more peaks)

We are now going to use maxima.p for the alignement. First we use the file ana0.par.

#########################################
### Really important parameters
## to be set to discover what the geometry conventions
## really are


orientation_codes=7  # 8 possibilities look at analyse_maxima_fourier.py
theta =  0.0
theta_offset = 0.0

alpha = 50.0
beta = 0.0
angular_step = +0.1

det_origin_X = 1242.06
det_origin_Y = 1258.6

dist =  299.987
pixel_size = 0.172
kappa = 0
omega = 90


# lmbda is not so important in the desperate phase. It just expand the size of the lattice
lmbda = 0.6968

#####################################################
### Other parameters that are generally small or that can
##  be refined later like r1,r2,r3
beam_tilt_angle = -0.05118

d1 = 0.04845
d2 = 0.1531

omega_offset =  0
phi = 0


##################################
# AA,BB,CC r1,r2,r3 : we will find them later 
############################"
## aAA, aBB, aCC... later


################################################
###### peaks positions are read from this file
maxima_file = "maxima.p"


###########################################################################
## we are going to visualise the spots distribution in reciprocal 3D space
##
transform_only  = 1

In this file we have put a rough knowledge of the experimental geometry and we have not set r1,r2,r3 yet, that must be aligned, and cell parameters neither. The first step is a visual check of the experimental geometry, which is triggered by setting, in ana0.par; transform_only to 1, so that q3d.txt will be generated by the following command:

tds2el_analyse_maxima_fourier_v1 ana0.par

and then visually inspected by

tds2el_rendering_v1 q3.txt
_images/calcite_0.png

The regular distribution of peaks is a sign that the geometry is correct. If you change ordering_codes ( the orientation of the detector) you may happen to find another choice of the orientation which works correctly, here the alternative would be 5. We do here the procedure with choice 7. Otherwise you can have slightly wrong parameters and still have a usuable prealignement. In this case an error of 10 or more pixels on the detector zero, would not be dramatic and you can use a rough idea of the geometry to bootstrap and refine later the uncertain parameters.

Given that the alignement is not so bad rerun with this this. input file . Fourier alignement is now performed, due to r1,r2,r3 set to zero*which triggers the alignement of the sample. Here how the input file looks like

#########################################
### Really important parameters
## to be set to discover what the geometry conventions
## really are

maxima_file = "maxima.p"

orientation_codes=7  # 8 possibilities look at analyse_maxima_fourier.py
theta =  0.0
theta_offset = 0.0

alpha = 50.0
beta = 0.0
angular_step = +0.1

det_origin_X = 1242.06
det_origin_Y = 1258.6

dist =  299.987
pixel_size = 0.172
kappa = 0
omega = 90

# lmbda is not so important in the desperate phase. It just expand the size of the lattice
lmbda = 0.6968

#####################################################
### Other parameters that are generally small or that can
##  be refined later like r1,r2,r3
beam_tilt_angle = -0.05118

d1 = 0.04845
d2 = 0.1531

omega_offset =  0
phi = 0

##################################
# AA,BB,CC r1,r2,r3 : settint r1,r2,r3 to zero
# triggers Fourier estimation of  r1,r2,r3,AA,BB,CC
# and also aAA,aBB,aCC
############################"
## 
r1 =  0
r2 = 0
r3 =  0
AA = 4.98
BB = 4.98
CC = 17.0

## this gives some constraint for the cell angles
aAA = 90
aBB = 90
aCC = 120



################################################
###### peaks positions are read from this file
maxima_file = "maxima.p"





The program stops producing the following estimations

#############
## ESTIMATION by FOURIER
AA=4.899338e+00
BB=4.899338e+00
CC=1.703103e+01
aAA=90.0            ##
aBB=90.0            ##  set by hand
aCC=120.0           ##
r3=-6.352717e+01
r2=6.700895e+01
r1=-1.735652e+02
############################

which is displayed on the screen. In reality several estimations are produced, ranked from the worst to the best according to the distance of the found AA,BB,CC from the values given in input. Copy (choose one if several are available) and paste a set of guessed parameters in the input file so that you get a new file for the next step. You can copy and past this estimation at the end of the ana1.par file, set now transform_only=1 and generate again q3d.txt to see how well it is aligned, remember the definition of GIACOVAZZI where astar is aligned with X axis.

In this run aAA,aBB,aCC where more or less constrained to their initial values and we have set them back to their nominal value. You can free one of them by setting its initial value to None. In all cases the search rejects combination of axis whose 3x3 determinant of their versors is smaller than 0.3(hard-coded parameter). The variable variations is an empty dictionary to signify that the program has to stop after producing the first guess of r1,r2,r3 and axes based on Fourier transform of the peaks distribution.

You can now run a fine retuning of the parameters by resetting transform_only to zero and adding variation informations for the choosed to-be-optimised variables (see the main documentation for explanations), we use the file ana2.par, obtained by adding first the initial guess for r1,r2,r3,AA,BB,CC and eventually aAA,aBB,aCC, plus the following intervals for variations at the end of the file, we have fixed there aAA, aBB, aCC by hand :: Note that in this optimisation above we have used dipendent AA and BB. You can set them to vary identically using this mechanism: you use an OrderedDict which keeps the ordering and, here below, BB which is set to -1 is interpreted to vary with the -1 in relative position, which is AA

variations =  collections.OrderedDict([
    ["r1", minimiser.Variable(0.0,-2,2)            ]    ,
    ["r2", minimiser.Variable(0,-6,10)             ]    ,
    ["r3", minimiser.Variable(0.0,-2,2)            ]     ,
    ["AA", minimiser.Variable(0.0,-0.1,0.1)        ]     ,
    ["BB",  -1   ]  ,
    ["CC", minimiser.Variable(0.0,-0.2,0.2)        ]      ,
  ]
)
## a tolerance. It controls when the fit stops
##
ftol = 1.0e-7

The relative position is taken on the ordered list of independent variables. Therefore for a cubic crystal you should set the variation like this

["AA", minimiser.Variable(0.0,-0.1,0.1)        ]     ,
["BB",  -1   ]  ,
["CC",  -1        ]      ,

because BB is not going in the list of independent variables, then -1 for CC means AA. In some cases it may be useful to add more variables to this dictionary, so that they will be optimised. In our particular case we know already well the geometry and we need to align just the sample.

There is also a more precise retuning of the geometry that we can do at Cij fitting time that we will see later in the documentation.

You can now run ana2.par. Now the optimisation is done with simplex descend, which is slow but roboust. You get refined parameters at the end, and at each step they are printed on screen. Reuse them for the next step, of course the first thing to do is to check if they look good, with transform_only=1. You know already how to do, you get this:

_images/calcite_1.png

At the end of the run a new maxima file tagged_maxima.p is written. Now you can select or discard spots by setting review=True and file_maxima=”tagged_maxima.p” in extraction.par and rerun tds2el_extraction. You generate newpeaks.p

_images/selected.png

But it is not necessary for now, we will work with the whole set of peaks. Now we are finally ready to go to the tds2el core. We start from the file ana3.par which contains the good geometry. On this basis we construct tds2el.par : which is formed adding few parameters to ana3.par (which already contains structural refined parameters from alignement process):

### To show a list of harvested spots and stop, uncomment line below
## show_hkls = True
##### manual selection, uncomment  line below
### collectSpots  = [[-5, 3, 1], [-5, 3, 4], [-5, 4, 0], [-4, 2, -6], [-4, 2, 0], [-4, 2, 6]]
#  to select spots from peaks file
file_maxima = "newpeaks.p"
### OPTIONS FOR EASY EDITING
DO_HARVESTING = 0
TODO = DO_HARVESTING
if(TODO==DO_HARVESTING):
    threshold = 100000.0
    DeltaQ_max = 0.15
    DeltaQ_min = 0.03
    load_from_file= None
    save_to_file  =  "collection4c.h5"
    DO_ELASTIC=1

and structuring in such a way that you can add new tasks and select them through TODO variable. The threshold parameters is used to filter out all those images with pixel values greater than threshold, while computing for each spot its center which is the weighted average of all points with value exceeding threshold. You can change this behaviour setting variable badradius to greater than zero, for example

badradius=300

to discard not the whole image but a disc of radius 300 around the bad pixels. If badradius==0 then all the image is discarded if only it contains one hot pixel. The harvested file collection4c.h5 will contain these baricenter informations, while the Qs are calculated based on geometry. Then we create the main input file tds2el.inp

conf_file             =  "tds2el.par"  
filter_file           =  "data/filter.edf"
images_prefix         =  "data/Calcite_170K_1_0"
ascii_generators      =  "2(110),3(z),1~"
normalisation_file    = "normalisation.txt"

We use the normalisation.txt that we have created with this script where we rely on the timestamps to interpolate in the esrf-machine file of the intensities. NOTE that we had to rely in the script on the files time stamps because all the metadata of one single scan files are all identical copy-and-paste. This might not be a problem though because for this Calcite example we are going to retrieve the relative elastic constants which depend mainly on the TDS shape around each peak, each peak being fitted with an ad-hic factor, and that a sweep through an hkl spot is relatively fast. Multiple sweep across the same hkl, we repeat it, are fitted separately with two independent factors.

When you run the harvesting, you can gain some time using mpi for parallel processing

mpirun -n 12 tds2el_v1 tds2el.inp

Then you have a lot of message flowing on the screen, you can check at least than in the lines like the one below, the hkl for the maximum of the image , at least the maximum with big value, are at close-to-integer values

the maximum is  1.15984e+06 [ 0.58050209  0.60439211 -0.23556644]  and hkl are  [ 0.99992228  3.01368737 -4.00607014]   fname  data/Calcite_170K_1_0113.cbf  i,j  (array([2000]), array([248]))

Now we redo a harvesting but a centered one this time, and we write it on collection_centered.h5 using the centering informations contained in collection4c.h5 to shift the spots blobs

DO_CENTERED_HARVESTING = 1
TODO = DO_CENTERED_HARVESTING
if(TODO==DO_CENTERED_HARVESTING):
    centering_file ="collection4c.h5"
    threshold = 100000.0
    DeltaQ_max = 0.15
    DeltaQ_min = 0.03
    load_from_file= None
    save_to_file  =  "collection_centered.h5"
    DO_ELASTIC=1

and rerun the harvesting obtaining collection_centered.h5.

Now do a dry run

DRY_RUN=2
if TODO == DRY_RUN:
    DO_ELASTIC=1
    load_from_file= "collection_centered.h5"
    save_to_file  =  None
    fitSpots = True
    file_maxima = "tagged_maxima.p"
    REPLICA_COLLECTION = 0

Getting the tensor structure

 GET ELASTIC SYMMETRY
1.000000 A11  | 1.000000 A11  | 1.000000 A13  |               | 1.000000 A15  |
              | -2.000000 A66 |               |               |               |
---------------------------------------------------------------------------------------------
1.000000 A11  | 1.000000 A11  | 1.000000 A13  |               | -1.000000 A15 |
-2.000000 A66 |               |               |               |               |
---------------------------------------------------------------------------------------------
1.000000 A13  | 1.000000 A13  | 1.000000 A33  |               |               |
              |               |               |               |               |
---------------------------------------------------------------------------------------------
              |               |               | 1.000000 A44  |               | -1.000000 A15
              |               |               |               |               |
---------------------------------------------------------------------------------------------
1.000000 A15  | -1.000000 A15 |               |               | 1.000000 A44  |
              |               |               |               |               |
---------------------------------------------------------------------------------------------
              |               |               | -1.000000 A15 |               | 1.000000 A66
              |               |               |               |               |
[(0, 0), (1, 1), (2, 2), (3, 3), (4, 4), (5, 5), (0, 1), (0, 2), (0, 3), (0, 4), (0, 5), (1, 2), (1, 3), (1, 4), (1, 5), (2, 3), (2, 4), (2, 5), (3, 4), (3, 5), (4, 5)]
 THE ELASTIC TENSOR VARIABLES ARE
['A11', 'A33', 'A44', 'A66', 'A13', 'A15']
 NOW STOPPING, TO DO TRHE FIT YOU NEED TO SET CGUESS
 processo  0  aspetta

And finally the last step

FIT=3
TODO = FIT
if TODO == FIT :
    DO_ELASTIC=1
    load_from_file= "collection_centered.h5"
    save_to_file  =  None
    fitSpots = True
    file_maxima = "tagged_maxima.p"

    CGUESS={"A11":  [13.7,0 , 400],
            "A13":  [2.37825229066,-400 ,400 ],
            "A15":  [-2.37825229066,-400 ,400 ],
            "A33":  [7.97 ,-400 ,400 ],
            "A44":  [3.42,-400 ,400 ],
            "A66":  [2.37825229066,-400 ,400 ],
    }

    doMinimisation=1
    visSim=1

    testSpot_beta = 0.06
    testSpot_Niter = 80
    testSpot_N = 128

that you launch NOT in parallel because it goes on GPU (or openCL CPU)

tds2el_v1 tds2el.inp

When we run the here provided tds2el.par input file (the three steps DO_HARVESTING, DO_CENTERED_HARVEST, DOFIT) we get at the last iteration

Cd  [  8.77377685  11.94596129   9.13998372  12.02890254  10.69581529
   9.06863885   6.83834762  12.24906892  13.60318685  12.41514349
  10.19839079  11.20609415   9.91145907  12.37931566   9.81276988
   9.96252389  10.57292079  13.71866018  12.35521767  13.19460209
  13.83431925  15.67267406  16.7445851   15.73214582  10.06162759
  11.62850902  13.98696931  18.20378866  16.82188007  14.67623455
  13.81360169  14.00076727   6.74455121  15.54030662  14.77989534
  10.8232562   16.64759301  16.53157886  21.34324891  10.36225485
  11.76831842  23.73767826  15.47190807  10.04734147  13.80688889
  16.05622385   8.74746071  12.50448148  14.79274612  11.55283805
   8.38569546  12.60616975  13.40084559  13.71701158   8.06683765
   9.62859946   9.59526193  12.72029506   8.74297481  13.20600079
  12.80713714  11.62122354  12.69188326  11.18902653   9.30684766
  10.60141131  10.32288793  10.1140334   11.45190703   7.86751104]
+------+---------------+---------------+---------------+---------------+----------------+-----+-----+
| A11  |      A33      |      A44      |      A66      |      A13      |      A15       |  r1 |  r2 |
+------+---------------+---------------+---------------+---------------+----------------+-----+-----+
| 15.6 | 8.57846854485 | 3.60973484091 | 4.84774053083 | 5.48408442006 | -2.11935386044 | 0.0 | 0.0 |
| 15.6 |     -400.0    |     -400.0    |     -400.0    |     -400.0    |     -400.0     | 0.0 | 0.0 |
| 15.6 |     400.0     |     400.0     |     400.0     |     400.0     |     400.0      | 0.0 | 0.0 |
+------+---------------+---------------+---------------+---------------+----------------+-----+-----+
+-----+-----+-----+-----+-------+-------+------+-------+
|  r3 |  AA |  BB |  CC | omega | alpha | beta | kappa |
+-----+-----+-----+-----+-------+-------+------+-------+
| 0.0 | 0.0 | 0.0 | 0.0 |  0.0  |  0.0  | 0.0  |  0.0  |
| 0.0 | 0.0 | 0.0 | 0.0 |  0.0  |  0.0  | 0.0  |  0.0  |
| 0.0 | 0.0 | 0.0 | 0.0 |  0.0  |  0.0  | 0.0  |  0.0  |
+-----+-----+-----+-----+-------+-------+------+-------+
+-----+-----+--------------+--------------+
|  d1 |  d2 | det_origin_X | det_origin_Y |
+-----+-----+--------------+--------------+
| 0.0 | 0.0 |     0.0      |     0.0      |
| 0.0 | 0.0 |     0.0      |     0.0      |
| 0.0 | 0.0 |     0.0      |     0.0      |
+-----+-----+--------------+--------------+
 error  617011795.024

The optimised parameters are shown in the table, with the current value of each placed above its interval. Here we have just optimised the elastic constants, but if we open the interval for the other we can simultaneaously perform a geometry optimisation. However here we have proceeded by recentering on Bragg peaks. A residual misalignement of the geometry could consist in a small rotation of the diffusions blobs around their Bragg peak and affect weakly the elastic constants. We will study later in the documentation a case where we do finer alignement, by opening the other intervals.

If you wonder what Cd is, it is the factor of the background values, a constant background for each spots. they are constrained variables of the optimisation, all as the per-spot factors which are invisible to the user, but multiply each calculated spot.

If you have set visSim=1 you can get at the end of the fit, if you are patient, back-interpolated volumes for each calculated and fitted spots. They are written in simtest.h5. To skim a compared view through selected slices of each spot

python skimma.py simetest.h5

using this script skimma.py.