Ticket #4513: RAWB.FOR;8

File RAWB.FOR;8, 13.6 KB (added by Peter Parker, 9 years ago)

Backward Scattering Angles Routine

Line 
1      PROGRAM DIFF
2      CHARACTER*80 runid
3      character runo*4,run1*5,fname*80,fname1*80
4      PARAMETER(NM=5,ND=200)
5
6      REAL T(10000),CT(10000),err(10000),SUM1(9),SUM2(9),SUM3(9)
7      REAL COUT(10000),CTHIN(10000),CTHICK(10000),NFR(9),sum4(9)
8      REAL ERROUT(10000),ERRTHIN(10000),ERRTHICK(10000)
9      REAL OUT(2000),THIN(2000),THICK(2000)
10      REAL MOUT(10000),MTHIN(10000),MTHICK(10000)
11      real xout(10000),yout(10000),eout(10000)
12      INTEGER IARRAY(100),ERRCODE,LOUT,NCHAN,NDET,IFSN,NOS
13      INTEGER rpb1(32),ICT(10000),ICTA(9,10000),MICTA(9,10000),NI(9)
14
15      write(6,*) ' first and last run numbers?'
16      read(5,*) irun1,irun2
17*      write(6,*) ' first and last spectrum numbers?'
18      ispectmin=3
19      ispectmax=134
20*      read(5,*) ispectmin,ispectmax
21*      write(6,*) ' for single difference type 1'
22*      write(6,*) ' for double difference type 2'
23*      read(5,*) iopt
24      iopt=2
25      beta=0.28
26*      write(6,*) ' enter file name for output'
27*      read(5,'(a)') fname
28      fname='tempb.dat'
29      fname1=fname
30*      write(6,'(80a)')fname1
31      open(unit=4,file=fname1,status='new')
32
33      do ISPECT=ispectmin,ispectmax    ! start of do loop over detectors
34      DO I=1,9
35      DO J=1,10000
36       ICTA(I,J)=0
37       MICTA(I,J)=0
38      END DO
39      END DO
40      NGOOD_FRAMES=0   
41      sum1=0.0
42      sum2=0.0
43
44      write(6,*) ' Spectrum; ',ispect
45      IF(IOPT.EQ.1) THEN
46*       WRITE(6,*) ' SINGLE DIFFERENCE'
47      ELSE IF(IOPT.EQ.2) THEN
48*       WRITE(6,*) ' DOUBLE DIFFERENCE'
49*       WRITE(6,*) ' BETA=',BETA
50      ELSE IF(IOPT.EQ.3) THEN
51       WRITE(6,*) ' THICK DIFFERENCE'
52      END IF
53
54
55      DO 101 IRUN=IRUN1,IRUN2
56* Define name of data file.
57*      write(6,*) ' irun=',irun
58      if(irun.lt.10000) then
59       write(runo(1:4),'(I4.4)')irun
60       runid='evs_data:evs0'//runo//'.raw'
61      else
62       write(run1(1:5),'(I5.5)')irun
63       runid='evs_data:evs'//run1//'.raw'
64      end if
65*      write(6,5) runid
66    5 Format(' ',' data file is;',50a)
67
68* Get relevant run information.
69      CALL GETPARI(runid,'RPB',rpb1,32,length_out,errcode)
70      ngood_frames =ngood_frames+rpb1(10)
71      CALL GETPARI(runid,'NPER',NPD,100,LOUT,ERRCODE)
72      CALL GETPARI(RUNID,'NSP1',IARRAY,100,LOUT,ERRCODE)
73      NDET=IARRAY(1)
74
75      CALL GETPARR(RUNID,'TCB1',T,10000,LOUT,ERRCODE)
76      CALL GETPARI(RUNID,'NTC1',IARRAY,100,LOUT,ERRCODE)
77      NCHAN=IARRAY(1)
78      NMON=MON(RUNID) 
79      DT=T(3)-T(2)
80      DO I=1,NCHAN
81       T(I)=T(I)-DT
82      END DO
83
84      DO IP=1,NPD
85       IS=ISPECT+(IP-1)*(NDET+1)
86       CALL GETDAT(RUNID,IS,1,ICT,10000,ERRCODE)
87       DO IT=2,NCHAN
88        ICTA(IP,IT)=ICTA(IP,IT)+ICT(IT)
89       END DO
90       IS=NMON+(IP-1)*(NDET+1)
91       CALL GETDAT(RUNID,IS,1,ICT,10000,ERRCODE)
92       DO IT=2,NCHAN
93        MICTA(IP,IT)=MICTA(IP,IT)+ICT(IT)
94       END DO
95      END DO
96
97       DO I=2,NCHAN,50
98        IF(I.LT.500.0) THEN
99*        write(6,*) i, i,icta(3,i),icta(4,i),icta(5,i),icta(6,i)
100        END IF
101       END DO
102
103
104  101 CONTINUE
105
106*      write(6,*) ' Total number of good frames=',ngood_frames
107*      WRITE(6,*) ' NUMBER OF PERIODS=',NPD
108*      WRITE(6,*) ' NUMBER OF SPECTRA=',NDET
109*      WRITE(6,*) ' Spectrum number=',ISPECT
110*      WRITE(6,*) ' Number of monitor spectrum=',nmon
111* DETERMINE WHICH SPECTRA IS OUT,THIN,THICK
112      IOUT=1
113      ITHIN=2
114      ITHICK=3
115
116       DO IP=1,NPD
117        DO IT=2,NCHAN
118         IF(T(IT).GT.300.0.AND.T(IT).LT.400.0) THEN
119          SUM1(IP)=SUM1(IP)+ICTA(IP,IT)
120         ELSE IF(T(IT).GT.550.0.AND.T(IT).LT.600.0) THEN
121          SUM2(IP)=SUM2(IP)+ICTA(IP,IT)
122         END IF
123        END DO
124*        write(6,*) ' sum1=',sum1(ip),' sum2=',sum2(ip)
125       END DO
126       
127       DO IP=1,NPD
128        NI(IP)=IP
129        IF(SUM2(IP).NE.0.0) THEN
130         SUM1(IP)=SUM1(IP)/SUM2(IP)
131         SUM4(IP)=SUM1(IP)
132        END IF
133       END DO
134
135*       do i=1,6
136*       write(6,*)i, ni(i),sum1(i)
137*       end do
138*       write(6,*)'2'
139
140       CALL ORDER(SUM1,NI,NPD)
141
142*       do i=1,6
143*       write(6,*)i, ni(i),sum1(i)
144*       end do
145
146       CALL WRITESP(NI,ISPECT,NDET,NPD) ! Write, out,thin,thick spectra.
147
148*      WRITE(6,*) ' 1'
149* Calculate total foil out, thin, thick and convert to Cts/micro-sec       
150      IF(NPD.EQ.2) THEN
151       DO I=2,NCHAN
152        COUT(I)=FLOAT(ICTA(NI(2),I))/(T(I)-T(I-1))
153        ERROUT(I)=SQRT(FLOAT(ICTA(NI(2),I)))/(T(I)-T(I-1))
154        CTHIN(I)=FLOAT(ICTA(NI(1),I))/(T(I)-T(I-1))
155        ERRTHIN(I)=SQRT(FLOAT(ICTA(NI(1),I)))/(T(I)-T(I-1))
156       END DO
157       SUM3(IOUT)=SUM2(NI(2))
158       SUM3(ITHIN)=SUM2(NI(1))
159       SUM3(ITHICK)=0.0
160      ELSE IF(NPD.EQ.3) THEN
161       DO I=2,NCHAN
162        COUT(I)=FLOAT(ICTA(NI(3),I))/(T(I)-T(I-1))
163        ERROUT(I)=SQRT(FLOAT(ICTA(NI(3),I)))/(T(I)-T(I-1))
164        CTHIN(I)=FLOAT(ICTA(NI(2),I))/(T(I)-T(I-1))
165        ERRTHIN(I)=SQRT(FLOAT(ICTA(NI(2),I)))/(T(I)-T(I-1))
166        CTHICK(I)=FLOAT(ICTA(NI(1),I))/(T(I)-T(I-1))
167        ERRTHICK(I)=SQRT(FLOAT(ICTA(NI(1),I)))/(T(I)-T(I-1))
168       END DO
169      SUM3(IOUT)=SUM2(NI(3))
170      SUM3(ITHIN)=SUM2(NI(2))
171      SUM3(ITHICK)=SUM2(NI(1))
172
173      ELSE IF(NPD.EQ.6) THEN
174       DO I=2,NCHAN
175        COUT(I)=FLOAT(ICTA(NI(5),I)+ICTA(NI(6),I))
176     $/(T(I)-T(I-1))
177        ERROUT(I)=SQRT(FLOAT(ICTA(NI(5),I)+ICTA(NI(6),I)))
178     $/(T(I)-T(I-1))
179        CTHIN(I)=FLOAT(ICTA(NI(3),I)+ICTA(NI(4),I)
180     $)/(T(I)-T(I-1))
181        ERRTHIN(I)=SQRT(FLOAT(ICTA(NI(3),I)+ICTA(NI(4),I)
182     $))/(T(I)-T(I-1))
183        CTHICK(I)=FLOAT(ICTA(NI(1),I)+ICTA(NI(2),I)
184     $)/(T(I)-T(I-1))
185        ERRTHICK(I)=SQRT(FLOAT(ICTA(NI(1),I)+ICTA(NI(2),I)
186     $))/(T(I)-T(I-1))
187       END DO
188      SUM3(IOUT)=SUM2(NI(5))+SUM2(NI(6))
189      SUM3(ITHIN)=SUM2(NI(3))+SUM2(NI(4))
190      SUM3(ITHICK)=SUM2(NI(1))+SUM2(NI(2))
191
192
193      ELSE IF(NPD.EQ.9) THEN
194       DO I=2,NCHAN
195        COUT(I)=FLOAT(ICTA(NI(7),I)+ICTA(NI(8),I)+
196     $ICTA(NI(9),I))/(T(I)-T(I-1))
197        ERROUT(I)=SQRT(FLOAT(ICTA(NI(7),I)+ICTA(NI(8),I)
198     $+ICTA(NI(9),I)))/(T(I)-T(I-1))
199        CTHIN(I)=FLOAT(ICTA(NI(4),I)+ICTA(NI(5),I)
200     $+ICTA(NI(6),I))/(T(I)-T(I-1))
201        ERRTHIN(I)=SQRT(FLOAT(ICTA(NI(4),I)+ICTA(NI(5),I)
202     $+ICTA(NI(6),I)))/(T(I)-T(I-1))
203        CTHICK(I)=FLOAT(ICTA(NI(1),I)+ICTA(NI(2),I)
204     $+ICTA(NI(3),I))/(T(I)-T(I-1))
205        ERRTHICK(I)=SQRT(FLOAT(ICTA(NI(1),I)+ICTA(NI(2),I)
206     $+ICTA(NI(3),I)))/(T(I)-T(I-1))
207       END DO
208      SUM3(IOUT)=SUM2(NI(7))+SUM2(NI(8))+SUM2(NI(9))
209      SUM3(ITHIN)=SUM2(NI(4))+SUM2(NI(5))+SUM2(NI(6))
210      SUM3(ITHICK)=SUM2(NI(1))+SUM2(NI(2))+SUM2(NI(3))
211      END IF
212
213       
214* Calculate total foil out, thin, thick monitor cts         
215      IF(NPD.EQ.2) THEN
216       DO I=2,NCHAN
217        MOUT(I)=FLOAT(MICTA(NI(2),I))/(T(I)-T(I-1))
218        MTHIN(I)=FLOAT(MICTA(NI(1),I))/(T(I)-T(I-1))
219       END DO
220      ELSE IF(NPD.EQ.3) THEN
221       DO I=2,NCHAN
222        MOUT(I)=FLOAT(MICTA(NI(3),I))/(T(I)-T(I-1))
223        MTHIN(I)=FLOAT(MICTA(NI(2),I))/(T(I)-T(I-1))
224        MTHICK(I)=FLOAT(MICTA(NI(1),I))/(T(I)-T(I-1))
225       END DO
226      ELSE IF(NPD.EQ.6) THEN
227       DO I=2,NCHAN
228        MOUT(I)=FLOAT(MICTA(NI(5),I)+MICTA(NI(6),I)
229     $)/(T(I)-T(I-1))
230        MTHIN(I)=FLOAT(MICTA(NI(3),I)+MICTA(NI(4),I)
231     $)/(T(I)-T(I-1))
232        MTHICK(I)=FLOAT(MICTA(NI(1),I)+MICTA(NI(2),I)
233     $)/(T(I)-T(I-1))
234       END DO
235      ELSE IF(NPD.EQ.9) THEN
236       DO I=2,NCHAN
237        MOUT(I)=FLOAT(MICTA(NI(7),I)+MICTA(NI(8),I)
238     $   +MICTA(NI(9),I))/(T(I)-T(I-1))
239        MTHIN(I)=FLOAT(MICTA(NI(4),I)+MICTA(NI(5),I)
240     $   +MICTA(NI(6),I))/(T(I)-T(I-1))
241        MTHICK(I)=FLOAT(MICTA(NI(1),I)+MICTA(NI(2),I)
242     $   +MICTA(NI(3),I))/(T(I)-T(I-1))
243       END DO
244      END IF     
245
246
247* Calculate number of frames in each period.
248      SUM3T=0.0
249      DO I=1,3
250       SUM3T=SUM3T+SUM3(I)
251      END DO
252      DO I=1,3
253       if(sum3t.ne.0.0)NFR(I)=NGOOD_FRAMES*SUM3(I)/SUM3T
254       IF(NFR(I).EQ.0.0) THEN
255         write(6,*) ' nframes=0 in period',i
256         nfr(i)=1
257       end if
258      END DO
259     
260
261* Start of routines for dead-time correction
262* CONVERT TO COUNTS/FRAME
263*      DO I=2,NCHAN
264*       COUT(I)=COUT(I)/NFR(IOUT)
265*       ERROUT(I)=ERROUT(I)/NFR(IOUT)
266*       CTHIN(I)=CTHIN(I)/NFR(ITHIN)
267*       ERRTHIN(I)=ERRTHIN(I)/NFR(ITHIN)
268*       IF(NPD.NE.2) THEN
269*        CTHICK(I)=CTHICK(I)/NFR(ITHICK)
270*        ERRTHICK(I)=ERRTHICK(I)/NFR(ITHICK)
271*       END IF
272*      END DO
273*
274* Correct for dead times in detector
275*      TAU=1.0
276*      DO I=2,NCHAN
277*       COUT(I)=COUT(I)/(1.0-TAU*COUT(I))
278*       ERROUT(I)=ERROUT(I)/(1.0-TAU*ERROUT(I))
279*       CTHIN(I)=CTHIN(I)/(1.0-TAU*CTHIN(I))
280*       ERRTHIN(I)=ERRTHIN(I)/(1.0-TAU*ERRTHIN(I))
281*       IF(NPD.NE.2) THEN
282*        ERRTHICK(I)=ERRTHICK(I)/(1.0-TAU*ERRTHICK(I))
283*       END IF
284*      END DO
285*
286* CONVERT BACK TO TOTAL COUNTS.
287*     DO I=2,NCHAN
288*       COUT(I)=COUT(I)*NFR(IOUT)
289*       ERROUT(I)=ERROUT(I)*NFR(IOUT)
290*       CTHIN(I)=CTHIN(I)*NFR(ITHIN)
291*       ERRTHIN(I)=ERRTHIN(I)*NFR(ITHIN)
292*       IF(NPD.NE.2) THEN
293*        CTHICK(I)=CTHICK(I)*NFR(ITHICK)
294*        ERRTHICK(I)=ERRTHICK(I)*NFR(ITHICK)
295*       END IF
296*      END DO
297
298
299
300* Calculate monitor cts between 600 and 700 usec.       
301      RMOUT=0.0
302      RMTHIN=0.0
303      RMTHICK=0.0
304      DO I=2,NCHAN
305       IF(T(I).GE.600.0.AND.T(I).LT.700.0) THEN
306        RMOUT=RMOUT+MOUT(I)
307        RMTHIN=RMTHIN+MTHIN(I)
308        RMTHICK=RMTHICK+MTHICK(I)
309       END IF
310      END DO
311      if(rmthick.eq.0.0)rmthick=1e-6
312* Normalise spectra to monitor counter.
313      const=1000.0 ! Constant to give suitable scale in tof spectra.
314      DO I=2,NCHAN
315       COUT(I)=COUT(I)*const/RMOUT
316       CTHIN(I)=CTHIN(I)*CONST/RMTHIN
317       CTHICK(I)=CTHICK(I)*CONST/RMTHICK
318       ERROUT(I)=ERROUT(I)*CONST/RMOUT
319       ERRTHIN(I)=ERRTHIN(I)*CONST/RMTHIN
320       ERRTHICK(I)=ERRTHICK(I)*CONST/RMTHICK
321      END DO
322
323
324* Calculate cts between 500 and 600 usec.       
325      SUMOUT=0.0
326      SUMTHIN=0.0
327      SUMTHICK=0.0
328      DO I=2,NCHAN
329       IF(T(I).GE.410.0.AND.T(I).LT.430.0) THEN
330        SUMOUT=SUMOUT+COUT(I)
331        SUMTHIN=SUMTHIN+CTHIN(I)
332        SUMTHICK=SUMTHICK+CTHICK(I)
333       END IF
334      END DO
335
336
337      if(sumthin.eq.0.0) then
338        write(6,*) ' no cts in thin foil spectra'
339        sumthin=1.0
340      end if
341      if(sumout.eq.0.0) then
342        write(6,*) ' no cts in foil out spectra'
343        sumout=1.0
344      end if
345      if(sumthick.eq.0.0) then
346        write(6,*) ' no cts in thick foil spectra'
347        sumthick=1.0
348      end if
349 
350* Normalise spectra to same area as foil out between 430 and 460 usec.
351      DO I=2,NCHAN
352       CTHIN(I)=CTHIN(I)*SUMOUT/SUMTHIN
353       ERRTHIN(I)=ERRTHIN(I)*SUMOUT/SUMTHIN
354       CTHICK(I)=CTHICK(I)*SUMOUT/SUMTHICK
355       ERRTHICK(I)=ERRTHICK(I)*SUMOUT/SUMTHICK
356      END DO
357
358      IT=0.0
359      IF(IOPT.EQ.1)THEN ! THIN DIFFERENCE
360       DT=T(3)-T(2)
361       DO I=2,NCHAN
362        IF(T(I).LT.600.0) THEN
363         IT=IT+1
364         XOUT(IT)=T(I)
365         YOUT(IT)=COUT(I)-CTHIN(I)
366         EOUT(IT)=SQRT(ERROUT(I)**2+ERRTHIN(I)**2)
367        END IF
368       END DO
369      ELSE IF(IOPT.EQ.2) THEN ! DOUBLE DIFFERENCE
370       DO I=2,NCHAN
371        IF(T(I).LT.600.0) THEN
372         IT=IT+1
373         XOUT(IT)=T(I)
374         YOUT(IT)=COUT(I)*(1-BETA)-CTHIN(I)+BETA*CTHICK(I)
375         EOUT(I)=SQRT((1-BETA)**2*ERROUT(I)**2+ERRTHIN(I)**2
376     $           +BETA**2*ERRTHICK(I)**2)
377        END IF
378       END DO
379      ELSE IF(IOPT.EQ.3) THEN ! THICK DIFFERENCE
380       DO I=2,NCHAN
381        IF(T(I).LT.600.0) THEN
382         IT=IT+1
383         XOUT(IT)=T(I)
384         YOUT(IT)=COUT(I)-CTHICK(I)
385         EOUT(IT)=SQRT(ERROUT(I)**2+ERRTHICK(I)**2)
386       END IF
387       END DO
388      END IF   
389
390*       DO I=2,NCHAN,50
391*        IF(T(I).LT.600.0) THEN
392*        write(6,*) t(i),cout(i),cthin(i),cthick(i)
393*        END IF
394*       END DO
395
396
397
398      LPTOUT=IT-1
399      write(4,*) ispect,lptout     
400      do i=1,lptout
401      write(4,*) xout(i),yout(i),eout(i)
402      end do
403
404      do i=1,lptout,20
405*      write(6,*) xout(i),yout(i),eout(i)
406      end do
407
408      end do ! end of do loop over detectors
409      write(4,*) irun1,irun2
410      close(4)
411      write(6,*) ' Data is in TEMPB.DAT'
412      END
413
414c Sort in order of increasing x values.
415      SUBROUTINE ORDER(X,Y,N)
416      REAL X(N)
417      INTEGER Y(N)
418
419      DO IS=1,N
420       IMIN=IS
421       XMIN=X(IS)
422       DO I=IS,N
423        IF(X(I).LE.XMIN) THEN
424        IMIN=I
425        XMIN=X(I)
426        END IF
427       END DO
428* Put smallest x value and corresponding y value in X(IS) and Y(IS).
429       XTEMP=X(IMIN)
430       YTEMP=Y(IMIN)
431       X(IMIN)=X(IS)
432       Y(IMIN)=Y(IS)
433       X(IS)=XTEMP
434       Y(IS)=YTEMP
435      END DO
436
437      END
438
439      SUBROUTINE WRITESP(NI,ISPECT,NDET,NPD)
440      INTEGER NI(9)
441       If(npd.eq.9) then
442       write(6,*) ' Foil out Spectra are:'
443       do i=7,9
444         IV=ISPECT+(ni(i)-1)*(NDET+1)
445         write(6,*) iv
446       end do
447
448       write(6,*) ' Thin Foil Spectra are:'
449       do i=4,6
450         IV=ISPECT+(ni(i)-1)*(NDET+1)
451         write(6,*) iv
452       end do
453
454       write(6,*) ' Thick foil Spectra are:'
455       do i=1,3
456         IV=ISPECT+(ni(i)-1)*(NDET+1)
457         write(6,*) iv
458       end do
459      end if
460
461       If(npd.eq.3) then
462       write(6,*) ' Foil out Spectrum is:'
463         IV=ISPECT+(ni(3)-1)*(NDET+1)
464         write(6,*) iv
465       write(6,*) ' Thin foil Spectrum is:'
466         IV=ISPECT+(ni(2)-1)*(NDET+1)
467         write(6,*) iv
468       write(6,*) ' Thick foil Spectrum is:'
469         IV=ISPECT+(ni(1)-1)*(NDET+1)
470         write(6,*) iv
471       end if
472
473       If(npd.eq.2) then
474       write(6,*) ' Foil out Spectrum is:'
475         IV=ISPECT+(ni(2)-1)*(NDET+1)
476         write(6,*) iv
477       write(6,*) ' Thin foil Spectrum is:'
478         IV=ISPECT+(ni(1)-1)*(NDET+1)
479         write(6,*) iv
480       end if
481       RETURN
482       END
483
484*  FUNCTION TO READ NUMBER OF MONITOR COUNTER
485      function mon(runid)
486      integer nmon*4,nrun*4
487      CHARACTER runid*80
488
489      CALL GETPARI(runid,'MDET',NMON,100,LOUT,ERRCODE)
490      MON=NMON
491      END