In [1]:
%%javascript
$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')


Goal: Slices through bhp_e

I learned that we can't go from bhp in time space to energy space because of the nonuniform binning, so I built a method to generate bhm_e on an event by event bases (methods > build_bhm_with_energy.ipynb), then plot as bhp_e (methods > build_plot_bhp_e.ipynb), and now I'm going to do some slices.

P. Schuster
University of Michigan
June 28 2018


In [1]:
import matplotlib.pyplot as plt
import matplotlib.colors
import numpy as np
import os
import scipy.io as sio
import sys
import time
import inspect
import pandas as pd
from tqdm import *

# Plot entire array
np.set_printoptions(threshold=np.nan)

In [2]:
%load_ext autoreload
%autoreload 2

In [3]:
sys.path.append('../scripts/')
import bicorr as bicorr
import bicorr_plot as bicorr_plot
import bicorr_e as bicorr_e
import bicorr_math as bicorr_math

Steps in the process:

  • Load the data for bhm_e
  • Produce bhp_e
  • Choose energies at which to take slices
  • Calculate slices, view them on bhp_e plot
  • Calculate energy averages

Load the data for bhm_e

Same data as methods > build_plot_bhp_e.ipynb.


In [4]:
bhm_e, e_bin_edges, note = bicorr_e.load_bhm_e('../analysis/Cf072115_to_Cf072215b/datap')

In [5]:
print(bhm_e.shape)
print(e_bin_edges.shape)
print(note)


(990, 1, 600, 600)
(601,)
Combined measurements from Cf072115 Cf072115b Cf072215a Cf072215b

In [6]:
det_df = bicorr.load_det_df()

In [7]:
dict_pair_to_index, dict_index_to_pair, dict_pair_to_angle = bicorr.build_dict_det_pair(det_df)

In [8]:
num_fissions = 2194651200.00

Produce bhp_e


In [11]:
num_fissions = 2194651200.00
bhp_e, norm_factor = bicorr_e.build_bhp_e(bhm_e,e_bin_edges,num_fissions=num_fissions,print_flag=True)


Creating bhp_e for...
pair_is =  all
energy bin width (MeV) =  0.025
length of pair_is =  990
norm_factor =  1357940430.0000002

In [12]:
bicorr_plot.bhp_e_plot(bhp_e, e_bin_edges, zoom_range=[0,6], title = "normalized bhp_e", show_flag = True)


C:\Users\pfschus\AppData\Local\Continuum\Anaconda3\lib\site-packages\matplotlib\cbook\deprecation.py:107: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.
  warnings.warn(message, mplDeprecation, stacklevel=1)
<Figure size 576x396 with 0 Axes>

Slices analysis

Choose energies at which to take slices


In [13]:
e_slices = list(np.arange(0.5,6,.5))
print(e_slices)


[0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5]

In [14]:
bicorr_plot.bhp_e_plot(bhp_e, e_bin_edges, zoom_range=[0,6], title='Plot of bhp_e', show_flag=False, clear_flag=False)

for e in e_slices:    
    plt.axvline(e,c='w') 
    
plt.show()


C:\Users\pfschus\AppData\Local\Continuum\Anaconda3\lib\site-packages\matplotlib\cbook\deprecation.py:107: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.
  warnings.warn(message, mplDeprecation, stacklevel=1)

Calculate slices


In [15]:
bhp_e_slice, slice_e_range = bicorr_e.slice_bhp_e(bhp_e, e_bin_edges, 2, 2.45, True)


Creating slice through bhp_e for energies from 2.0 to 2.475

In [16]:
bhp_e_slices,slice_e_ranges = bicorr_e.slices_bhp_e(bhp_e,e_bin_edges,e_slices,0.224)

In [17]:
len(e_bin_edges)


Out[17]:
601

In [18]:
bhp_e_slices.shape


Out[18]:
(11, 600)

In [19]:
slice_e_ranges


Out[19]:
array([[0.5  , 0.725],
       [1.   , 1.225],
       [1.5  , 1.725],
       [2.   , 2.225],
       [2.5  , 2.725],
       [3.   , 3.225],
       [3.5  , 3.725],
       [4.   , 4.225],
       [4.5  , 4.725],
       [5.   , 5.225],
       [5.5  , 5.725]])

Plot the slices


In [20]:
bicorr_plot.plot_bhp_e_slice(bhp_e_slice,e_bin_edges,slice_e_range,
                             show_flag=True)


<Figure size 576x396 with 0 Axes>

In [21]:
bicorr_plot.plot_bhp_e_slices(bhp_e_slices,e_bin_edges,slice_e_ranges)


Out[21]:
['0.50 to 0.73',
 '1.00 to 1.23',
 '1.50 to 1.73',
 '2.00 to 2.23',
 '2.50 to 2.73',
 '3.00 to 3.23',
 '3.50 to 3.73',
 '4.00 to 4.23',
 '4.50 to 4.73',
 '5.00 to 5.23',
 '5.50 to 5.73']

In [22]:
bicorr_plot.plot_bhp_e_slices(bhp_e_slices,e_bin_edges,slice_e_ranges,
                              show_flag = False)
plt.xlim([0,6])
plt.ylim([0,0.01])
plt.show()


Calculate energy averages

Need to work in absolute counts, not normalized counts.


In [23]:
bhp_e, norm_factor = bicorr_e.build_bhp_e(bhm_e,e_bin_edges,print_flag=True)


Creating bhp_e for...
pair_is =  all
energy bin width (MeV) =  0.025
length of pair_is =  990
norm_factor =  1

In [24]:
bhp_e_slices,slice_e_ranges = bicorr_e.slices_bhp_e(bhp_e,e_bin_edges,e_slices,0.224)
print(slice_e_ranges)


[[0.5   0.725]
 [1.    1.225]
 [1.5   1.725]
 [2.    2.225]
 [2.5   2.725]
 [3.    3.225]
 [3.5   3.725]
 [4.    4.225]
 [4.5   4.725]
 [5.    5.225]
 [5.5   5.725]]

In [25]:
E_min = 0.75; E_max = 4;

Eave, Eave_err, Ej = bicorr_e.calc_Eave_slices(bhp_e_slices,e_slices,e_bin_edges,E_min,E_max)

Plot it.


In [36]:
bicorr_plot.plot_Eave_vs_Ej(Eave, Eave_err, Ej)


Over what range was this taken?


In [37]:
bicorr_plot.plot_bhp_e_slices(bhp_e_slices,e_bin_edges,slice_e_ranges,
                              E_min = E_min, E_max = E_max,
                              show_flag = False)
plt.xlim([0,6])
plt.ylim([0,0.01])
plt.show()


In 10-degree bins


In [9]:
# Set up data

In [10]:
angle_bin_edges = np.arange(8,190,10)
angle_bin_centers = bicorr_math.calc_centers(angle_bin_edges)
angle_indices = [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,17]

In [11]:
e_slices = list(np.arange(0.5,6,.5))
E_min = 0.75; E_max = 4;

In [12]:
# Allocate arrays

In [13]:
bhp_e = np.zeros((len(angle_bin_centers),len(e_bin_edges)-1,len(e_bin_edges)-1))
norm_factor = np.zeros(len(angle_bin_centers))

bhp_e_slices = np.zeros((len(angle_bin_centers),len(e_slices),len(e_bin_edges)-1))

Eave = np.zeros((len(angle_bin_centers),len(e_slices)))
Eave_err = np.zeros((len(angle_bin_centers),len(e_slices)))

In [30]:
# Do the calculations
for i in angle_indices: #range(len(angle_bin_centers)):
    angle_min = angle_bin_edges[i]
    angle_max = angle_bin_edges[i+1]
    
    pair_is = bicorr.generate_pair_is(det_df, angle_min, angle_max)
    bhp_e[i,:,:], norm_factor[i] = bicorr_e.build_bhp_e(bhm_e,e_bin_edges,pair_is=pair_is,num_fissions = num_fissions,print_flag=True)
    bhp_e_slices[i,:,:],slice_e_ranges = bicorr_e.slices_bhp_e(bhp_e[i,:,:],e_bin_edges,e_slices,0.224)
    Eave[i,:], Eave_err[i,:], _ = bicorr_e.calc_Eave_slices(bhp_e_slices[i,:,:],e_slices,e_bin_edges,E_min,E_max,norm_factor=norm_factor[i])


Creating bhp_e for...
pair_is =  [  0   8  44  87 129 170 210 249 287 332 360 368 395 429 462 494 525 555
 584 620 639 665 690 714 737 759 780 800 837 845 854 870 885 899 912 924
 935 953 954 962 969 975 980 984 987 989]
energy bin width (MeV) =  0.025
length of pair_is =  46
norm_factor =  63096222.000000015
Creating bhp_e for...
pair_is =  [  9  51  52 295 296 331 369 402 403 592 593 619 846 861 862 943 944 952]
energy bin width (MeV) =  0.025
length of pair_is =  18
norm_factor =  24689826.000000004
Creating bhp_e for...
pair_is =  [  1  10  17  45  53  88  93  94  95  96 130 136 137 138 171 177 178 179
 211 217 218 219 250 256 257 258 259 294 330 341 361 370 396 404 430 435
 436 437 438 463 469 470 471 495 501 502 503 526 532 533 534 556 562 563
 564 565 591 618 640 656 666 691 715 738 760 781 836 838 847 855 863 871
 876 877 878 879 886 892 893 894 900 906 907 908 913 919 920 921 925 931
 932 933 934 942 951 955 963 970 976 981 985 988]
energy bin width (MeV) =  0.025
length of pair_is =  102
norm_factor =  139909014.00000003
Creating bhp_e for...
pair_is =  [  2  18  46  54  60  61  89  97 131 135 139 172 176 180 212 216 220 255
 293 304 305 340 362 397 405 431 439 464 468 472 496 500 504 527 531 535
 561 590 641 667 692 716 739 761 839 856 864 872 880 887 891 895 901 905
 909 914 918 922 930 941 956 964 971 977 982 986]
energy bin width (MeV) =  0.025
length of pair_is =  66
norm_factor =  90529362.00000001
Creating bhp_e for...
pair_is =  [ 11  19  35  55  62  98 102 103 104 134 140 175 181 215 221 254 266 267
 268 292 303 329 339 359 371 377 386 406 440 467 473 499 505 530 536 560
 589 617 629 638 647 657 681 682 817 818 827 835 848 865 881 890 896 904
 910 917 923 929 940 950]
energy bin width (MeV) =  0.025
length of pair_is =  60
norm_factor =  82299420.00000001
Creating bhp_e for...
pair_is =  [  3  12  20  26  47  56  63  90  99 105 132 141 143 144 145 146 147 173
 174 182 186 187 188 214 226 227 228 229 230 253 265 291 302 328 338 350
 363 372 387 398 407 420 432 441 465 474 497 498 506 529 559 588 611 616
 637 642 648 658 668 672 683 693 705 706 717 740 798 799 809 816 826 834
 840 849 857 866 873 882 888 897 902 903 911 916 928 939 949 957 965 972
 978 983]
energy bin width (MeV) =  0.025
length of pair_is =  92
norm_factor =  126192444.00000003
Creating bhp_e for...
pair_is =  [  4  27  36  48  64  69  78  91 106 107 133 148 149 184 185 189 190 224
 225 263 264 301 314 323 349 358 364 378 388 399 411 421 433 453 466 583
 602 610 628 636 643 649 659 669 673 684 694 696 707 708 718 728 729 730
 731 776 777 778 779 790 796 797 808 815 825 833 841 858 874 889 958 966
 973 979]
energy bin width (MeV) =  0.025
length of pair_is =  74
norm_factor =  101502618.00000001
Creating bhp_e for...
pair_is =  [ 13  21  37  57  65  79 100 108 120 142 150 183 191 213 223 252 262 286
 290 300 322 327 337 357 373 379 408 412 422 442 444 454 475 528 558 574
 582 587 601 609 615 627 660 674 685 697 709 732 733 750 751 752 753 754
 755 756 757 758 774 775 789 795 807 814 832 850 867 883 898 915 927 938
 948]
energy bin width (MeV) =  0.025
length of pair_is =  73
norm_factor =  100130961.00000001
Creating bhp_e for...
pair_is =  [  5  14  22  28  49  58  66  70  92 101 109 111 151 222 251 261 277 289
 299 313 326 336 348 365 374 389 400 409 423 434 443 455 485 486 553 554
 557 581 586 608 614 635 644 650 670 675 686 695 698 710 711 719 720 734
 735 769 770 772 773 788 793 794 806 813 824 842 851 859 868 875 884 926
 937 947 959 967 974]
energy bin width (MeV) =  0.025
length of pair_is =  77
norm_factor =  105617589.00000001
Creating bhp_e for...
pair_is =  [  6  23  29  38  50  67  71  80 110 112 121 152 161 239 248 260 276 285
 298 312 321 335 347 356 366 380 390 401 413 424 445 456 476 487 516 517
 523 524 545 552 573 580 600 607 626 634 645 651 661 671 676 687 699 712
 721 736 741 742 748 749 768 771 787 792 805 812 823 831 843 860 960 968]
energy bin width (MeV) =  0.025
length of pair_is =  72
norm_factor =  98759304.00000001
Creating bhp_e for...
pair_is =  [ 15  24  39  59  68  81 122 162 201 209 247 284 288 297 320 325 334 355
 375 381 391 410 414 425 446 457 477 488 489 490 492 493 507 515 518 519
 520 521 522 544 546 547 549 550 551 572 579 585 599 606 613 625 633 652
 662 677 688 700 713 722 723 724 726 727 743 744 745 746 747 762 763 765
 766 767 786 791 804 811 822 830 852 869 936 946]
energy bin width (MeV) =  0.025
length of pair_is =  84
norm_factor =  115219188.00000003
Creating bhp_e for...
pair_is =  [  7  16  30  72  82 113 123 153 163 192 200 202 208 238 246 275 283 311
 319 324 346 367 376 392 415 426 447 458 459 460 461 478 491 508 514 543
 548 571 575 576 577 578 598 605 612 632 646 653 663 678 689 701 702 703
 704 725 764 782 783 784 785 803 810 821 829 844 853 945 961]
energy bin width (MeV) =  0.025
length of pair_is =  69
norm_factor =  94644333.00000001
Creating bhp_e for...
pair_is =  [ 25  31  40  73 114 124 154 160 164 169 193 199 203 207 231 237 240 245
 274 282 310 333 345 354 382 427 448 479 484 509 513 537 542 570 604 624
 664 679 802 828]
energy bin width (MeV) =  0.025
length of pair_is =  40
norm_factor =  54866280.00000001
Creating bhp_e for...
pair_is =  [ 41  83 125 128 165 167 168 204 205 206 241 242 244 278 281 318 353 383
 393 416 428 449 452 480 482 483 510 511 512 538 539 541 566 569 597 603
 623 631 654 680 801 820]
energy bin width (MeV) =  0.025
length of pair_is =  42
norm_factor =  57609594.000000015
Creating bhp_e for...
pair_is =  [ 32  42  74  84  85  86 115 119 126 127 155 159 166 194 198 232 236 243
 269 273 279 280 309 315 316 317 344 352 384 394 417 418 419 450 451 481
 540 567 568 594 595 596 622 630 655 819]
energy bin width (MeV) =  0.025
length of pair_is =  46
norm_factor =  63096222.000000015
Creating bhp_e for...
pair_is =  [ 33  43  75  77 116 118 156 158 195 197 233 235 270 272 306 308 343 351
 385 621]
energy bin width (MeV) =  0.025
length of pair_is =  20
norm_factor =  27433140.000000004
Creating bhp_e for...
pair_is =  [ 34  76 117 157 196 234 271 307 342]
energy bin width (MeV) =  0.025
length of pair_is =  9
norm_factor =  12344913.000000002

In [31]:
vmin = np.min(bhp_e[np.nonzero(bhp_e)])
vmax = np.max(bhp_e)

In [32]:
Eave_min = np.min(Eave[np.nonzero(Eave)])
Eave_max = np.max(Eave)

In [36]:
# Make the plots
filenames_bhp_e = []
filenames_Eave  = []  
for i in angle_indices: #range(len(angle_bin_centers)):
    angle_min = angle_bin_edges[i]
    angle_max = angle_bin_edges[i+1]   
    
    title = '{} to {} degrees'.format(angle_min, angle_max)  
    filename_bhp_e = 'bhp_e_{}_{}_deg'.format(angle_min,angle_max); filenames_bhp_e.append(filename_bhp_e);
    bicorr_plot.bhp_e_plot(bhp_e[i,:,:], e_bin_edges, zoom_range = [0,6], 
                           vmin=vmin, vmax=vmax,
                           title=title, show_flag = True,
                           save_flag = True, save_filename = filename_bhp_e)        
    #bicorr_plot.plot_bhp_e_slices(bhp_e_slices[i,:,:],e_bin_edges,slice_e_ranges,
    #                               E_min = E_min, E_max = E_max, title=title,
    #                               save_filename = 'bhp_e_slices_{}_{}_degrees'.format(angle_min,angle_max))
    filename_Eave = 'Eave_{}_{}_degrees'.format(angle_min,angle_max); filenames_Eave.append(filename_Eave);
    bicorr_plot.plot_Eave_vs_Ej(Eave[i,:], Eave_err[i,:], e_slices, title=title,
                                y_range = [Eave_min,Eave_max],
                                save_flag = True, save_filename = filename_Eave)


<Figure size 576x396 with 0 Axes>
<Figure size 576x396 with 0 Axes>
<Figure size 576x396 with 0 Axes>
<Figure size 576x396 with 0 Axes>
<Figure size 576x396 with 0 Axes>
<Figure size 576x396 with 0 Axes>
<Figure size 576x396 with 0 Axes>
<Figure size 576x396 with 0 Axes>
<Figure size 576x396 with 0 Axes>
<Figure size 576x396 with 0 Axes>
<Figure size 576x396 with 0 Axes>
<Figure size 576x396 with 0 Axes>
<Figure size 576x396 with 0 Axes>
<Figure size 576x396 with 0 Axes>
<Figure size 576x396 with 0 Axes>
<Figure size 576x396 with 0 Axes>
<Figure size 576x396 with 0 Axes>

Save data to file


In [38]:
np.savez('datap/slices_analysis', 
          angle_bin_edges = angle_bin_edges,
          angle_bin_centers = angle_bin_centers,
          angle_indices = angle_indices,
          e_slices = e_slices,
          E_min = E_min, E_max = E_max,
          bhp_e = bhp_e, norm_factor=norm_factor,
          bhp_e_slices = bhp_e_slices,
          Eave=Eave, Eave_err = Eave_err
          )

Animate it all


In [39]:
import imageio

In [42]:
images_bhp_e = []
for filename in filenames_bhp_e:
    images_bhp_e.append(imageio.imread(os.path.join('fig',filename + '.png')))
imageio.mimsave('fig/animate_bhp_e.gif',images_bhp_e, fps=1)

In [45]:
images_Eave = []
for filename in filenames_Eave:
    images_Eave.append(imageio.imread(os.path.join('fig',filename + '.png')))
imageio.mimsave('fig/animate_Eave.gif',images_Eave, fps=1)