-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathanalysis.pgevo.waterbox.py
executable file
·152 lines (103 loc) · 7.51 KB
/
analysis.pgevo.waterbox.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
#!/usr/bin/env python
import image,numpy as np,auger,plot,dump
from scipy.ndimage.filters import gaussian_filter
###########################################################################################################
volume_offset=-149#waterbox sources
pgsrc_ct = image.image('stage2/data/source-139.mhd')
pgsrc_ct.toprojection(".x", [0,1,1,1])
pgsrc_ct_y = pgsrc_ct.imdata/pgsrc_ct.imdata.max()
pgsrc_rpct = image.image('stage2/data/source-144.mhd')
pgsrc_rpct.toprojection(".x", [0,1,1,1])
pgsrc_rpct_y = pgsrc_rpct.imdata/pgsrc_rpct.imdata.max()
pgsrc_ct_x = np.linspace(-149,149,150) #bincenters
pgsrc_ct_xhist = np.linspace(-150,150,151) #2mm voxels, endpoints
pgsrc_ct_x = pgsrc_ct_x+(volume_offset-pgsrc_ct_x[0]) #offset for pg source image
pgsrc_ct_xhist = pgsrc_ct_xhist+(volume_offset-pgsrc_ct_xhist[0]) #same
pgsrc_rpct_x = pgsrc_ct_x #is nu hetzelfde
pgsrc_rpct_xhist = pgsrc_ct_xhist
pgsrc_ctfo=auger.get_fop(pgsrc_ct_x,pgsrc_ct_y)
pgsrc_rpctfo=auger.get_fop(pgsrc_rpct_x,pgsrc_rpct_y)
###########################################################################################################
#emission in worldframe
smooth_param = 8.5 #IBA: 20 FWHM of some spread function = 2.35sigma of a gaussian, Priegnitz 2015
pgemis_ct = dump.thist2np_xy('fromcluster/run.8ZQJ/output.2399064/GammProdCount-Prod.root')
pgemis_rpct = dump.thist2np_xy('fromcluster/run.ahDK/output.2399464/GammProdCount-Prod.root')
pgemis_ct_x = np.array(pgemis_ct['histo'][0])
pgemis_ct_y = np.array(pgemis_ct['histo'][1])
pgemis_ct_y = pgemis_ct_y/pgemis_ct_y.max()
pgemis_rpct_x = np.array(pgemis_rpct['histo'][0])
pgemis_rpct_y = np.array(pgemis_rpct['histo'][1])
pgemis_rpct_y = pgemis_rpct_y/pgemis_rpct_y.max()
pgemis_ct_y_smooth = gaussian_filter(pgemis_ct_y, sigma=smooth_param)
pgemis_ct_y_smooth = pgemis_ct_y_smooth/pgemis_ct_y_smooth.max()
pgemis_rpct_y_smooth = gaussian_filter(pgemis_rpct_y, sigma=smooth_param)
pgemis_rpct_y_smooth = pgemis_rpct_y_smooth/pgemis_rpct_y_smooth.max()
pgemis_ctfo=auger.get_fop(pgemis_ct_x,pgemis_ct_y)
pgemis_rpctfo=auger.get_fop(pgemis_rpct_x,pgemis_rpct_y)
pgemis_smooth_ctfo=auger.get_fop(pgemis_ct_x,pgemis_ct_y_smooth)
pgemis_smooth_rpctfo=auger.get_fop(pgemis_rpct_x,pgemis_rpct_y_smooth)
###########################################################################################################
#detection in local frame, remove shift. coords already centered as zero, so can just add the shift
det_offset = -17.31
ipnl = auger.getctset(1e9,'fromcluster/run.8ZQJ','fromcluster/run.ahDK','ipnl-auger-tof-1.root')
pgipnl_ct_x = np.array(ipnl['ct']['x'])
pgipnl_ct_x = pgipnl_ct_x+det_offset
pgipnl_ct_y = np.array(ipnl['ct']['av'])
pgipnl_ct_y = pgipnl_ct_y-(pgipnl_ct_y[pgipnl_ct_y < np.percentile(pgipnl_ct_y, 25)].mean()) #remove floor
pgipnl_ct_y = pgipnl_ct_y/pgipnl_ct_y.max() #scale
pgipnl_rpct_x = np.array(ipnl['rpct']['x'])
pgipnl_rpct_x = pgipnl_rpct_x+det_offset
pgipnl_rpct_y = np.array(ipnl['rpct']['av'])
pgipnl_rpct_y = pgipnl_rpct_y-(pgipnl_rpct_y[pgipnl_rpct_y < np.percentile(pgipnl_rpct_y, 25)].mean()) #remove floor
pgipnl_rpct_y = pgipnl_rpct_y/pgipnl_rpct_y.max() #scale
pgipnl_ctfo = np.mean(ipnl['ct']['falloff'])+det_offset
pgipnl_rpctfo = np.mean(ipnl['rpct']['falloff'])+det_offset
iba = auger.getctset(1e9,'fromcluster/run.8ZQJ','fromcluster/run.ahDK','iba-auger-notof-3.root')
pgiba_ct_x = np.array(iba['ct']['x'])
pgiba_ct_x = pgiba_ct_x+det_offset
pgiba_ct_y = np.array(iba['ct']['av'])
pgiba_ct_y = pgiba_ct_y-(pgiba_ct_y[pgiba_ct_y < np.percentile(pgiba_ct_y, 25)].mean()) #remove floor
pgiba_ct_y = pgiba_ct_y/pgiba_ct_y.max() #scale
pgiba_rpct_x = np.array(iba['rpct']['x'])
pgiba_rpct_x = pgiba_rpct_x+det_offset
pgiba_rpct_y = np.array(iba['rpct']['av'])
pgiba_rpct_y = pgiba_rpct_y-(pgiba_rpct_y[pgiba_rpct_y < np.percentile(pgiba_rpct_y, 25)].mean()) #remove floor
pgiba_rpct_y = pgiba_rpct_y/pgiba_rpct_y.max() #scale
pgiba_ctfo = np.mean(iba['ct']['falloff'])+det_offset
pgiba_rpctfo = np.mean(iba['rpct']['falloff'])+det_offset
###########################################################################################################
f, (ax1,ax2) = plot.subplots(nrows=1, ncols=2, sharex=False, sharey=False)
#ax1.step(pgsrc_ct_x,pgsrc_ct_y, color='steelblue',lw=1., alpha=0.2, label='PGSRC FOP: '+str(pgsrc_ctfo), where='mid')
#ax1.step(pgsrc_rpct_x,pgsrc_rpct_y, color='indianred',lw=1., alpha=0.2, label='PGSRC RPCT FOP: '+str(pgsrc_rpctfo), where='mid')
ax1.step(pgemis_ct_x,pgemis_ct_y, color='steelblue',lw=1., alpha=0.6, label='Prod FOP:\n'+str(pgemis_ctfo)[:5], where='mid')
#ax1.step(pgemis_rpct_x,pgemis_rpct_y, color='indianred',lw=1., alpha=0.4, label='PROD RPCT FOP: '+str(pgemis_rpctfo), where='mid')
ax1.step(pgemis_ct_x,pgemis_ct_y_smooth, color='seagreen',lw=1., alpha=0.6, label='Prod Smooth FOP:\n'+str(pgemis_smooth_ctfo)[:5], where='mid')
#ax1.step(pgemis_rpct_x,pgemis_rpct_y_smooth, color='indianred',lw=1., alpha=0.6, label='PROD SMOOTH RPCT FOP: '+str(pgemis_smooth_rpctfo), where='mid')
pgipnl_ctfo = auger.get_fop(pgipnl_ct_x,pgipnl_ct_y,plot=True,ax=ax1,label='Detected Smooth ')
#ax1.step(pgipnl_ct_x,pgipnl_ct_y, color='indianred',lw=1., alpha=0.8, label='PGIPNL FOP:\n'+str(pgipnl_ctfo)[:5], where='mid')
#ax1.step(pgipnl_rpct_x,pgipnl_rpct_y, color='indianred',lw=1., alpha=0.8, label='PGIPNL RPCT FOP: '+str(pgipnl_rpctfo), where='mid')
# CT WITHOUT SMOOTHING
ax2.step(pgemis_ct_x,pgemis_ct_y, color='steelblue',lw=1., alpha=0.6, label='Prod FOP:\n'+str(pgemis_ctfo)[:5], where='mid')
#ax1.step(pgemis_rpct_x,pgemis_rpct_y, color='indianred',lw=1., alpha=0.4, label='PROD RPCT FOP: '+str(pgemis_rpctfo), where='mid')
ax2.step(pgemis_ct_x,pgemis_ct_y_smooth, color='seagreen',lw=1., alpha=0.6, label='Prod Smooth FOP:\n'+str(pgemis_smooth_ctfo)[:5], where='mid')
#ax1.step(pgemis_rpct_x,pgemis_rpct_y_smooth, color='indianred',lw=1., alpha=0.6, label='PROD SMOOTH RPCT FOP: '+str(pgemis_smooth_rpctfo), where='mid')
pgipnl_ctfo = auger.get_fop(pgipnl_ct_x,pgipnl_ct_y,plot=True,ax=ax2,smooth=False,label='Detected ')
#ax2.step(pgipnl_ct_x,pgipnl_ct_y, color='indianred',lw=1., alpha=0., label='PGIPNL FOP:\n'+str(pgipnl_ctfo)[:5], where='mid')
#ax1.step(pgipnl_rpct_x,pgipnl_rpct_y, color='indianred',lw=1., alpha=0.8, label='PGIPNL RPCT FOP: '+str(pgipnl_rpctfo), where='mid')
########################################### rpct:
##ax2.step(pgsrc_ct_x,pgsrc_ct_y, color='steelblue',lw=1., alpha=0.2, label='PGSRC FOP: '+str(pgsrc_ctfo), where='mid')
##ax2.step(pgsrc_rpct_x,pgsrc_rpct_y, color='indianred',lw=1., alpha=0.2, label='PGSRC RPCT FOP: '+str(pgsrc_rpctfo), where='mid')
#ax2.step(pgemis_ct_x,pgemis_ct_y, color='steelblue',lw=1., alpha=0.4, label='PROD FOP: '+str(pgemis_ctfo), where='mid')
##ax2.step(pgemis_rpct_x,pgemis_rpct_y, color='indianred',lw=1., alpha=0.4, label='PROD RPCT FOP: '+str(pgemis_rpctfo), where='mid')
#ax2.step(pgemis_ct_x,pgemis_ct_y_smooth, color='steelblue',lw=1., alpha=0.6, label='PROD SMOOTH FOP: '+str(pgemis_smooth_ctfo), where='mid')
##ax2.step(pgemis_rpct_x,pgemis_rpct_y_smooth, color='indianred',lw=1., alpha=0.6, label='PROD SMOOTH RPCT FOP: '+str(pgemis_smooth_rpctfo), where='mid')
#ax2.step(pgiba_ct_x,pgiba_ct_y, color='steelblue',lw=1., alpha=0.8, label='PGIBA FOP: '+str(pgiba_ctfo), where='mid')
##ax2.step(pgiba_rpct_x,pgiba_rpct_y, color='indianred',lw=1., alpha=0.8, label='PGIBA RPCT FOP: '+str(pgiba_rpctfo), where='mid')
ax1.set_xlim(-50,50)
ax2.set_xlim(-50,50)
ax1.legend(frameon = False,loc='upper right')
ax2.legend(frameon = False,loc='upper right')
plot.texax(ax1)
plot.texax(ax2)
f.savefig('pgevo-waterbox.pdf', bbox_inches='tight')
plot.close('all')