1 jones 1.1 subroutine sf_lookup_init(infile,protonflag)
2
3 ! Read in spectral function. First line has to have # of Pm,Em bins.
4 ! Subsequent lines (number of lines = numPm*numEm) are of the form:
5 !
6 ! <Pmvalue> <Emvalue> <SF_PROB> <dPm> <dEm>
7 !
8 ! where Pmvalue/Emvalue are the center of the Pm/Em bins, dPm/dEm are widths,
9 ! and <SF_PROB> is the integrated spectral function within that bin,
10 ! (SF_PROB = SF * dE * P**2 * dP * 4 * pi), and is the probability within
11 ! the Em/Pm bin, IF the spectral function is normalized to 1.
12 !
13 implicit none
14
15 include 'sf_lookup.inc'
16 include 'target.inc'
17 include 'constants.inc'
18
19 real*8 tmpPm,tmpEm,tmpsfp,tmpsfn,tmpdEm,tmpdPm
20 integer*4 iPm,iEm
21 character infile*80
22 jones 1.1 logical protonflag
23
24 open(unit=55,file=infile,status='old')
25
26 read (55,*) numPm , numEm
27
28 do iPm = 1 , numPm
29 do iEm = 1 , numEm
30 read (55,*) tmpPm,tmpEm,tmpsfp,tmpsfn,tmpdPm,tmpdEm
31
32 ! Choose proton or neutron spectral function.
33 if (protonflag) then
34 sfval(iEm,iPm)=tmpsfp
35 else
36 sfval(iEm,iPm)=tmpsfn
37 endif
38
39 ! Check that Pm value is consistent with others in same bin
40 if (iEm.eq.1) then
41 Pmval(iPm)=tmpPm
42 dPm(iPm)=tmpdPm
43 jones 1.1 else
44 if (abs(tmpPm-Pmval(iPm)).gt.0.01) then
45 write(6,*) 'Bad Pm value, line number is approx',(numEm*(iPm-1)+iEm)
46 write(6,*) 'iPm,oldPm,newPm=',iPm,Pmval(iPm),tmpPm
47 endif
48 endif
49
50 ! Check that Em value is consistent with others in same bin
51 if (iPm.eq.1) then
52 Emval(iEm)=tmpEm
53 dEm(iEm)=tmpdEm
54 else
55 if (abs(tmpEm-Emval(iEm)).gt.0.01) then
56 write(6,*) 'Bad Em value, line number is approx',(numEm*(iPm-1)+iEm)
57 write(6,*) 'iEm,oldEm,newEm=',iEm,Emval(iEm),tmpEm
58 endif
59 endif
60 enddo
61 enddo
62
63 sftotnorm = 0.0
64 jones 1.1 do iPm = 1 , numPm
65 sfnorm(iPm) = 0.0
66 do iEm = 1 , numEm
67 sfnorm(iPm) = sfnorm(iPm) + sfval(iEm,iPm)
68 sftotnorm = sftotnorm + sfval(iEm,iPm)
69 enddo
70 ! write(6,*) 'Pm=',Pmval(iPm),' has normalization=',sfnorm(iPm)
71 enddo
72 write(6,*) 'Uncorrected S.F. has normalization=',sftotnorm
73
74 do iPm = 1 , numPm
75 do iEm = 1 , numEm
76 sfval(iEm,iPm) = sfval(iEm,iPm) / sftotnorm
77 enddo
78 enddo
79
80 return
81 end
82
83
84 !----------------------------------------------------------------------
85 jones 1.1
86 subroutine sf_lookup(Em,Pm,SF)
87
88 ! Get value of spectral function <SF_PROB> at Em, Pm.
89 !
90 ! Where <SF_PROB> is the integrated spectral function within that bin,
91 ! (SF_PROB = SF * dE * P**2 * dP * 4 * pi), and is the probability within
92 ! the Em/Pm bin, once the spectral function has been normalized to 1.
93 !
94
95 implicit none
96
97 include 'sf_lookup.inc'
98
99 integer*4 ind,iPm,iEm
100 real*8 Em,Pm
101 real*8 Em1,Em2,sf1,sf2,logsf,w1,w2
102 real*8 SF
103
104 !find nearest Pm value
105 if (Pm.ge.Pmval(numPm)) then
106 jones 1.1 iPm=numPm-1
107 w1=0
108 w2=1
109 else if (Pm.le.Pmval(1)) then
110 iPm=1
111 w1=1
112 w2=0
113 else
114
115 ind=1
116 do while(Pm.gt.Pmval(ind))
117 ind = ind + 1
118 enddo
119 iPm=ind
120 w2=(Pm-Pmval(iPm))/(Pmval(iPm+1)-Pmval(iPm))
121 w1=(Pmval(iPm+1)-Pm)/(Pmval(iPm+1)-Pmval(iPm))
122 if (abs(w1*Pmval(iPm)+w2*Pmval(iPm+1)-Pm).gt.0.0001) then
123 write(6,*) 'w1,w2,Pm,Pmval(iPm)=',w1,w2,Pm,Pmval(iPm)
124 stop
125 endif
126 ! write(6,*) 'Pm,iPm,Pm,w1,w2=',Pm,iPm,Pmval(iPm),w1,w2
127 jones 1.1
128 ! do ind=1,numPm
129 ! if (Pm.gt.Pmval(ind)) then
130 ! iPm=ind
131 ! w2=(Pm-Pmval(iPm))/(Pmval(iPm+1)-Pmval(iPm))
132 ! w1=(Pmval(iPm+1)-Pm)/(Pmval(iPm+1)-Pmval(iPm))
133 ! if (abs(w1*Pmval(iPm)+w2*Pmval(iPm+1)-Pm).gt.0.0001) then
134 ! write(6,*) 'w1,w2,Pm,Pmval(iPm)=',w1,w2,Pm,Pmval(iPm)
135 ! stop
136 ! endif
137 ! write(6,*) 'Pm,iPm,Pm,w1,w2=',Pm,iPm,Pmval(iPm),w1,w2
138 ! endif
139 ! enddo
140
141 endif
142 if (abs(w1+w2-1).gt.0.0001) then
143 write(6,*) 'iPm,Pm,w1+w2=',iPm,Pm,w1+w2
144 stop
145 endif
146
147 if (Em.le.Emval(1)) then !linear extrapolation of LOG(f(Em)).
148 jones 1.1 Em1=Emval(1)
149 Em2=Emval(2)
150 sf1=w1*sfval(1,iPm)+w2*sfval(1,iPm+1)
151 sf2=w1*sfval(2,iPm)+w2*sfval(2,iPm+1)
152
153 else if (Em.gt.Emval(numEm)) then !linear extrapolation of LOG(f(Em))
154 Em1=Emval(numEm-1)
155 Em2=Emval(numEm)
156 sf1=w1*sfval(numEm-1,iPm)+w2*sfval(numEm-1,iPm+1)
157 sf2=w1*sfval(numEm,iPm)+w2*sfval(numEm,iPm+1)
158
159 else
160 do iEm=1,numEm-1
161 if (Em.ge.Emval(iEm) .and. Em.lt.Emval(iEm+1)) then
162 Em1=Emval(iEm)
163 Em2=Emval(iEm+1)
164 sf1=w1*sfval(iEm,iPm)+w2*sfval(iEm,iPm+1)
165 sf2=w1*sfval(iEm+1,iPm)+w2*sfval(iEm+1,iPm+1)
166 ! write(6,*) 'iEm,w1,w2=',iEm,w1,w2
167 endif
168 enddo
169 jones 1.1 endif
170
171 logsf=(sf1 + (Em-Em1)*(sf2-sf1)/(Em2-Em1))
172 SF=logsf
173
174 if (SF.lt.1.e-20) SF=0
175
176 return
177 end
178
179
180 !----------------------------------------------------------------------
181
182 subroutine generate_em(Pm,Em)
183
184 ! Generate a missing energy for an event with a given Pm.
185 !
186 ! Have to get spectral function for several Em bins, at the given value
187 ! of Pm. Then, select Em according to the energy distribution at fixed Pm.
188 !
189
190 jones 1.1 implicit none
191
192 include 'sf_lookup.inc'
193
194 integer*4 ind,iEm
195 real*8 Pm,Em !input Pm value, generated Em value.
196 real*8 x(nEmmax),y(nEmmax)
197 real*8 ynorm
198 real*8 ranprob
199
200 real*8 grnd
201
202 ! Determine Em distribution at fixed Pm (integrated SF up to Em).
203 ynorm = 0.
204 x(1) = Emval(1)
205 call sf_lookup(x(1),Pm,y(1))
206
207 do iEm = 2 , numEm
208 x(iEm) = Emval(iEm)
209 call sf_lookup(x(iEm),Pm,y(iEm))
210 y(iEm) = y(iEm) + y(iEm-1)
211 jones 1.1 enddo
212
213 ! Normalize the Em distribution.
214 do iEm = 1 , numEm
215 y(iEm) = y(iEm) / y(numEm)
216 enddo
217
218 ! Generate Em.
219 ranprob = grnd()
220 ind = 1
221 do while (ranprob.gt.y(ind))
222 ind = ind + 1
223 enddo
224 Em = Emval(ind) + dEm(ind)*(grnd()-0.5)
225
226 ! CHECK NORMALIZATION, HAVE CODE SAVE EM/PM DISTRIBUTIONS TO CHECK.
227 ! MAKE SURE EM IS NOT BELOW EM_MIN.
228
229 ! write(99,'(2f8.2)') Pm,Em
230
231 return
232 jones 1.1 end
|
Return to
sf_lookup.f
CVS log
Up to [HallC] / Poltar