/[PAMELA software]/DarthVader/ToFLevel2/src/tofl2com.for
ViewVC logotype

Diff of /DarthVader/ToFLevel2/src/tofl2com.for

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1.1.1 by mocchiut, Sat Jun 17 12:14:56 2006 UTC revision 1.13 by mocchiut, Mon Nov 23 09:50:50 2009 UTC
# Line 1  Line 1 
1  *****************************************************************************  
2    C******************************************************************************
3    C
4    C  08-12-06 WM: adc_c-bug :  The raw ADc value was multiplied with cos(theta)
5    C  and AFTER that there was an if statement "if tof32(right,i,iadc) < 4095"
6    C
7    C  jan-07 GF: ADCflags(4,12) inserted to flag artificial ADC values
8    C  jan-07 WM: artificial ADC values created using attenuation calibration
9    C  jan-07 WM: modified xtofpos flag "101". xtofpos must be inside physical
10    C             dimension of the paddle +/- 10 cm
11    C  jan-07 WM: if xtofpos=101 then this paddle is not used for beta
12    C             calculation
13    C  jan-07 WM: the definition for a "hit" is changed: Now we must have a
14    C             valid TDC signal on both sides
15    C  jan-07 WM: flag for PMTs #10 and #35 added, TDC=819 due to bit-shift
16    C  jan-07 WM: bug removed: in some cases tdc_tw was calculated due to a
17    C             leftover "xhelp" value
18    C  apr-07 WM: attenuation fit curve is now a double exponential fit
19    C             conversion from raw ADC to pC using calibration function
20    C             variables xtr_tof and ytr_tof inserted (filled with default)
21    C  jan-08 WM: Major Update: Time Walk correction introduced
22    C             Additionalyl we use the information from the "check_charge"
23    C             function to fill artificial ADC values and make small corrections
24    C             to the k1-parameter (for Z>2)
25    C  feb-08 WM: Calculation of beta(13) changed: First a mean beta is calculated,
26    C             then in a second step we check the residuals of the single
27    C             measurements, reject if > 10 sigma, calculate chi2 and "quality"
28    C             beta is taken as good if chi2<20 and quality>10
29    C  mar-08 WM: Call to "newbeta" changed, now a flag tells the function if the
30    C             call comes from "tofl2com" or form "toftrack"
31    C  mar-08 WM: Bug found in dEdx if check_charge>1
32    C  oct-08 WM: Calculation of zenith angle debugged, sometimes strange values
33    C             were possible
34    C  nov-09 WM: the dEdx part ("adctof_c") moved to the new dEdx routine from Napoli
35    C******************************************************************************
36    
37        INTEGER FUNCTION TOFL2COM()        INTEGER FUNCTION TOFL2COM()
38  c      c    
39        IMPLICIT NONE        IMPLICIT NONE
# Line 13  C     Line 48  C    
48        LOGICAL check        LOGICAL check
49        REAL secure        REAL secure
50    
51        INTEGER j        INTEGER j,hitvec(6)
       REAL xhelp_a,xhelp_t  
52    
53        REAL dx,dy,dr,ds        REAL dx,dy,dr,ds
54        REAL yhelp,xdummy,xkorr0,xhelp,xhelp1,xhelp2        REAL yhelp,yhelp1,yhelp2,xhelp,xhelp1,xhelp2
55        REAL c1,c2,sw,sxw,w_i        REAL c1,c2
56        INTEGER icount  
57    C      REAL sw,sxw,w_i
58    C      INTEGER icount
59    C      REAL beta_mean
60    
61        INTEGER tof11_j,tof21_j,tof31_j        INTEGER tof11_j,tof21_j,tof31_j
62        INTEGER tof12_j,tof22_j,tof32_j        INTEGER tof12_j,tof22_j,tof32_j
63    
       REAL beta_mean  
   
   
64  c     value for status of each PM-data  c     value for status of each PM-data
65  c     first index  : 1 = left, 2 = right  c     first index  : 1 = left, 2 = right
66  c     second index : 1... number of paddle  c     second index : 1... number of paddle
# Line 34  c     second index : 1... number of padd Line 68  c     second index : 1... number of padd
68        INTEGER tof21_event(2,2),tof22_event(2,2)        INTEGER tof21_event(2,2),tof22_event(2,2)
69        INTEGER tof31_event(2,3),tof32_event(2,3)        INTEGER tof31_event(2,3),tof32_event(2,3)
70    
71    
72          REAL y_coor_lin11c(8,2),x_coor_lin12c(6,2)
73          REAL x_coor_lin21c(2,2),y_coor_lin22c(2,2)
74          REAL y_coor_lin31c(3,2),x_coor_lin32c(3,2)
75    
76          DATA y_coor_lin11c(1,1),y_coor_lin11c(1,2) /-20.66,-2.497/
77          DATA y_coor_lin11c(2,1),y_coor_lin11c(2,2) /-9.10, -2.52/
78          DATA y_coor_lin11c(3,1),y_coor_lin11c(3,2) /-24.07,-2.12/
79          DATA y_coor_lin11c(4,1),y_coor_lin11c(4,2) /-13.40,-2.47/
80          DATA y_coor_lin11c(5,1),y_coor_lin11c(5,2) /-31.07,-2.32/
81          DATA y_coor_lin11c(6,1),y_coor_lin11c(6,2) /-21.69,-2.63/
82          DATA y_coor_lin11c(7,1),y_coor_lin11c(7,2) /-12.37,-2.65/
83          DATA y_coor_lin11c(8,1),y_coor_lin11c(8,2) /-10.81,-3.15/
84    
85          DATA x_coor_lin12c(1,1),x_coor_lin12c(1,2) /12.96, -2.65/
86          DATA x_coor_lin12c(2,1),x_coor_lin12c(2,2) /17.12,-2.44/
87          DATA x_coor_lin12c(3,1),x_coor_lin12c(3,2) /7.26, -1.98/
88          DATA x_coor_lin12c(4,1),x_coor_lin12c(4,2) /-22.52,-2.27/
89          DATA x_coor_lin12c(5,1),x_coor_lin12c(5,2) /-18.54,-2.28/
90          DATA x_coor_lin12c(6,1),x_coor_lin12c(6,2) /-7.67,-2.15/
91    
92          DATA x_coor_lin21c(1,1),x_coor_lin21c(1,2) /22.56,-1.56/
93          DATA x_coor_lin21c(2,1),x_coor_lin21c(2,2) /13.94,-1.56/
94    
95          DATA y_coor_lin22c(1,1),y_coor_lin22c(1,2) /-24.24,-2.23/
96          DATA y_coor_lin22c(2,1),y_coor_lin22c(2,2) /-45.99,-1.68/
97    
98          DATA y_coor_lin31c(1,1),y_coor_lin31c(1,2) /-22.99,-3.54/
99          DATA y_coor_lin31c(2,1),y_coor_lin31c(2,2) /-42.28,-4.10/
100          DATA y_coor_lin31c(3,1),y_coor_lin31c(3,2) /-41.29,-3.69/
101    
102          DATA x_coor_lin32c(1,1),x_coor_lin32c(1,2) /0.961, -3.22/
103          DATA x_coor_lin32c(2,1),x_coor_lin32c(2,2) /4.98,-3.48/
104          DATA x_coor_lin32c(3,1),x_coor_lin32c(3,2) /-22.08,-3.37/
105    
106                
107        REAL theta12,theta13,theta23        REAL theta13
108  C--   DATA ZTOF/53.74,53.04,23.94,23.44,-23.49,-24.34/ !Sergio 9.05.2006  C--   DATA ZTOF/53.74,53.04,23.94,23.44,-23.49,-24.34/ !Sergio 9.05.2006
109        REAL tofarm12        REAL tofarm12
110        PARAMETER (tofarm12 = 29.70) ! from 53.39 to 23.69        PARAMETER (tofarm12 = 29.70) ! from 53.39 to 23.69
111        REAL tofarm23                REAL tofarm23        
112        PARAMETER (tofarm23 = 47.61)  ! from 23.69 to -23.92        PARAMETER (tofarm23 = 47.61)  ! from 23.69 to -23.92
113        REAL tofarm13        REAL tofarm13
114        PARAMETER (tofarm13 = 77.31)  ! from 53.39 to -23.92        PARAMETER (tofarm13 = 77.31)  ! from 53.39 to -23.92
115          
116          REAL hepratio      
117    
118        INTEGER ihelp        INTEGER ihelp
119        REAL xkorr        REAL xkorr,btemp(12)
120    
121          REAL atten,pc_adc,check_charge,newbeta
122    
123          INTEGER IZ
124          REAL k1corrA1,k1corrB1,k1corrC1
125    
126    
127          INTEGER ifst
128          DATA  ifst /0/
129    
130  C---------------------------------------  C---------------------------------------
131  C      C    
# Line 56  C     Line 135  C    
135  C      C    
136  C     CALCULATE COMMON VARIABLES  C     CALCULATE COMMON VARIABLES
137  C      C    
138    C-------------------------------------------------------------------
139    
140  *******************************************************************          if (ifst.eq.0) then
141        icounter = icounter + 1  
142            ifst=1
143    
144  *     amplitude has to be 'secure' higher than pedestal for an adc event  C     amplitude has to be 'secure' higher than pedestal for an adc event
145        secure = 2.         secure = 2.
146    
147    C     ratio between helium and proton ca. 4
148           hepratio = 4.  !
149           offset = 1
150           slope = 2
151           left = 1
152           right = 2
153           none_ev = 0
154           none_find = 0
155           tdc_ev = 1
156           adc_ev = 1
157           itdc = 1
158           iadc = 2
159    
160    C--- These are the corrections to the k1-value for Z>2 particles
161           k1corrA1 = 0.
162           k1corrB1 = -5.0
163           k1corrC1=  8.0
164    
165    
166            ENDIF
167    C---------------------------------------------------------------------
168    
169          icounter = icounter + 1
170    
       offset = 1  
       slope = 2  
       left = 1  
       right = 2  
       none_ev = 0  
       none_find = 0  
       tdc_ev = 1  
       adc_ev = 1  
       itdc = 1  
       iadc = 2  
171    
172        do i=1,13        do i=1,13
173           betatof_a(i) = 100.          ! As in "troftrk.for"           betatof_a(i) = 100.          ! As in "troftrk.for"
174        enddo        enddo
175    
176          do i=1,6
177             hitvec(i) = -1
178          enddo
179    
180        do i=1,4        do i=1,4
181           do j=1,12           do j=1,12
182              adctof_c(i,j) = 1000.              adctof_c(i,j) = 1000.
# Line 92  C     Line 191  C    
191        enddo        enddo
192    
193    
194  c the calibration files are read in the main program from xxx_tofcalib.rz        do i=1,12
195             do j=1,4
196                tofmask(j,i) = 0
197             enddo
198          enddo
199    
200    
201    c gf adc falg:
202          do i=1,4
203             do j=1,12
204                adcflagtof(i,j) = 0
205             enddo
206          enddo
207    
208    c gf tdc falg:
209          do i=1,4
210             do j=1,12
211                tdcflagtof(i,j) = 0
212             enddo
213          enddo
214    
215    
216    C---  Fill xtr_tof and ytr_tof: positions from tracker at ToF layers
217    C---  since this is standalone ToF fill with default values
218          do j=1,6
219          xtr_tof(j) = 101.
220          ytr_tof(j) = 101.
221          enddo
222    
223    c the calibration files are read in the main program from xxx_tofcalib.rz
224    
225  c-------------------------get ToF data --------------------------------  c-------------------------get ToF data --------------------------------
226    
227  c     put the adc and tdc values from ntuple into tofxx(i,j,k) variables  c     put the adc and tdc values from ntuple into tofxx(i,j,k) variables
228    c     adc valueas are then pC
229    
230        do j=1,8        do j=1,8
231           tof11(1,j,2) = adc(ch11a(j),hb11a(j))           tof11(1,j,2) = pc_adc(adc(ch11a(j),hb11a(j)))
232           tof11(2,j,2) = adc(ch11b(j),hb11b(j))           tof11(2,j,2) = pc_adc(adc(ch11b(j),hb11b(j)))
233           tof11(1,j,1) = tdc(ch11a(j),hb11a(j))           tof11(1,j,1) = (tdc(ch11a(j),hb11a(j)))
234           tof11(2,j,1) = tdc(ch11b(j),hb11b(j))           tof11(2,j,1) = (tdc(ch11b(j),hb11b(j)))
235        enddo        enddo
236    
237    
238        do j=1,6        do j=1,6
239           tof12(1,j,2) = adc(ch12a(j),hb12a(j))           tof12(1,j,2) = pc_adc(adc(ch12a(j),hb12a(j)))
240           tof12(2,j,2) = adc(ch12b(j),hb12b(j))           tof12(2,j,2) = pc_adc(adc(ch12b(j),hb12b(j)))
241           tof12(1,j,1) = tdc(ch12a(j),hb12a(j))           tof12(1,j,1) = (tdc(ch12a(j),hb12a(j)))
242           tof12(2,j,1) = tdc(ch12b(j),hb12b(j))           tof12(2,j,1) = (tdc(ch12b(j),hb12b(j)))
243        enddo        enddo
244    
245        do j=1,2        do j=1,2
246           tof21(1,j,2) = adc(ch21a(j),hb21a(j))           tof21(1,j,2) = pc_adc(adc(ch21a(j),hb21a(j)))
247           tof21(2,j,2) = adc(ch21b(j),hb21b(j))           tof21(2,j,2) = pc_adc(adc(ch21b(j),hb21b(j)))
248           tof21(1,j,1) = tdc(ch21a(j),hb21a(j))           tof21(1,j,1) = (tdc(ch21a(j),hb21a(j)))
249           tof21(2,j,1) = tdc(ch21b(j),hb21b(j))           tof21(2,j,1) = (tdc(ch21b(j),hb21b(j)))
250        enddo        enddo
251    
252        do j=1,2        do j=1,2
253           tof22(1,j,2) = adc(ch22a(j),hb22a(j))           tof22(1,j,2) = pc_adc(adc(ch22a(j),hb22a(j)))
254           tof22(2,j,2) = adc(ch22b(j),hb22b(j))           tof22(2,j,2) = pc_adc(adc(ch22b(j),hb22b(j)))
255           tof22(1,j,1) = tdc(ch22a(j),hb22a(j))           tof22(1,j,1) = (tdc(ch22a(j),hb22a(j)))
256           tof22(2,j,1) = tdc(ch22b(j),hb22b(j))           tof22(2,j,1) = (tdc(ch22b(j),hb22b(j)))
257        enddo        enddo
258    
259        do j=1,3        do j=1,3
260           tof31(1,j,2) = adc(ch31a(j),hb31a(j))           tof31(1,j,2) = pc_adc(adc(ch31a(j),hb31a(j)))
261           tof31(2,j,2) = adc(ch31b(j),hb31b(j))           tof31(2,j,2) = pc_adc(adc(ch31b(j),hb31b(j)))
262           tof31(1,j,1) = tdc(ch31a(j),hb31a(j))           tof31(1,j,1) = (tdc(ch31a(j),hb31a(j)))
263           tof31(2,j,1) = tdc(ch31b(j),hb31b(j))           tof31(2,j,1) = (tdc(ch31b(j),hb31b(j)))
264        enddo        enddo
265    
266        do j=1,3        do j=1,3
267           tof32(1,j,2) = adc(ch32a(j),hb32a(j))           tof32(1,j,2) = pc_adc(adc(ch32a(j),hb32a(j)))
268           tof32(2,j,2) = adc(ch32b(j),hb32b(j))           tof32(2,j,2) = pc_adc(adc(ch32b(j),hb32b(j)))
269           tof32(1,j,1) = tdc(ch32a(j),hb32a(j))           tof32(1,j,1) = (tdc(ch32a(j),hb32a(j)))
270           tof32(2,j,1) = tdc(ch32b(j),hb32b(j))           tof32(2,j,1) = (tdc(ch32b(j),hb32b(j)))
271        enddo        enddo
272    
273  C----------------------------------------------------------------------  C----------------------------------------------------------------------
# Line 188  C--------------------------------------- Line 315  C---------------------------------------
315           if (abs(tof32(2,i,iadc)).gt.10000.) tof32(2,i,iadc)= 10000.           if (abs(tof32(2,i,iadc)).gt.10000.) tof32(2,i,iadc)= 10000.
316        ENDDO        ENDDO
317    
318    C----------------------------------------------------------------------
319    C------------------  set ADC & TDC flag = 0    ------------------------
320    C----------------------------------------------------------------------
321    
322          do j=1,8
323          if (adc(ch11a(j),hb11a(j)).LT.4096)adcflagtof(ch11a(j),hb11a(j))=0
324          if (adc(ch11b(j),hb11b(j)).LT.4096)adcflagtof(ch11b(j),hb11b(j))=0
325          if (tdc(ch11a(j),hb11a(j)).LT.4096)tdcflagtof(ch11a(j),hb11a(j))=0
326          if (tdc(ch11b(j),hb11b(j)).LT.4096)tdcflagtof(ch11b(j),hb11b(j))=0
327          enddo
328          do j=1,6
329          if (adc(ch12a(j),hb12a(j)).LT.4096)adcflagtof(ch12a(j),hb12a(j))=0
330          if (adc(ch12b(j),hb12b(j)).LT.4096)adcflagtof(ch12b(j),hb12b(j))=0
331          if (tdc(ch12a(j),hb12a(j)).LT.4096)tdcflagtof(ch12a(j),hb12a(j))=0
332          if (tdc(ch12b(j),hb12b(j)).LT.4096)tdcflagtof(ch12b(j),hb12b(j))=0
333          enddo
334          do j=1,2
335          if (adc(ch21a(j),hb21a(j)).LT.4096)adcflagtof(ch21a(j),hb21a(j))=0
336          if (adc(ch21b(j),hb21b(j)).LT.4096)adcflagtof(ch21b(j),hb21b(j))=0
337          if (tdc(ch21a(j),hb21a(j)).LT.4096)tdcflagtof(ch21a(j),hb21a(j))=0
338          if (tdc(ch21b(j),hb21b(j)).LT.4096)tdcflagtof(ch21b(j),hb21b(j))=0
339          enddo
340          do j=1,2
341          if (adc(ch22a(j),hb22a(j)).LT.4096)adcflagtof(ch22a(j),hb22a(j))=0
342          if (adc(ch22b(j),hb22b(j)).LT.4096)adcflagtof(ch22b(j),hb22b(j))=0
343          if (tdc(ch22a(j),hb22a(j)).LT.4096)tdcflagtof(ch22a(j),hb22a(j))=0
344          if (tdc(ch22b(j),hb22b(j)).LT.4096)tdcflagtof(ch22b(j),hb22b(j))=0
345          enddo
346          do j=1,3
347          if (adc(ch31a(j),hb31a(j)).LT.4096)adcflagtof(ch31a(j),hb31a(j))=0
348          if (adc(ch31b(j),hb31b(j)).LT.4096)adcflagtof(ch31b(j),hb31b(j))=0
349          if (tdc(ch31a(j),hb31a(j)).LT.4096)tdcflagtof(ch31a(j),hb31a(j))=0
350          if (tdc(ch31b(j),hb31b(j)).LT.4096)tdcflagtof(ch31b(j),hb31b(j))=0
351          enddo
352          do j=1,3
353          if (adc(ch32a(j),hb32a(j)).LT.4096)adcflagtof(ch32a(j),hb32a(j))=0
354          if (adc(ch32b(j),hb32b(j)).LT.4096)adcflagtof(ch32b(j),hb32b(j))=0
355          if (tdc(ch32a(j),hb32a(j)).LT.4096)tdcflagtof(ch32a(j),hb32a(j))=0
356          if (tdc(ch32b(j),hb32b(j)).LT.4096)tdcflagtof(ch32b(j),hb32b(j))=0
357          enddo
358    
359  C----------------------------------------------------------------  C----------------------------------------------------------------
360  C------------Check Paddles for hits  -----------------------  C---------- Check PMTs 10 and 35 for strange TDC values----------
361    C----------------------------------------------------------------
362    
363    C---- S116A TDC=819
364           if (tof11(1,6,1).EQ.819) then
365                 tof11(1,6,1) = 4095
366                 tdcflagtof(ch11a(6),hb11a(6))=2
367           endif
368      
369    C---- S222B TDC=819
370           if (tof22(2,2,1).EQ.819) then
371                 tof22(2,2,1) = 4095
372                 tdcflagtof(ch22b(2),hb22b(2))=2
373           endif
374    
375    C----------------------------------------------------------------
376    C------------   Check Paddles for hits    -----------------------
377    C------  a "hit" means TDC values<4095 on both sides ------------
378  C----------------------------------------------------------------  C----------------------------------------------------------------
379    
380  C     upper tof  S11  C     upper tof  S11
# Line 199  C     upper tof  S11 Line 384  C     upper tof  S11
384              tof11_event(j,i) = none_ev              tof11_event(j,i) = none_ev
385              IF ((tof11(j,i,itdc).LT.2000).AND.(tof11(j,i,itdc).GT.100))              IF ((tof11(j,i,itdc).LT.2000).AND.(tof11(j,i,itdc).GT.100))
386       +           tof11_event(j,i) = tof11_event(j,i) + tdc_ev       +           tof11_event(j,i) = tof11_event(j,i) + tdc_ev
             IF ((tof11(j,i,iadc).GT.secure).AND.  
      +           (tof11(j,i,iadc).LT.4095))  
      +           tof11_event(j,i) = tof11_event(j,i) + adc_ev  
387           ENDDO           ENDDO
388        ENDDO        ENDDO
389    
# Line 232  C     upper tof  S12 Line 414  C     upper tof  S12
414              tof12_event(j,i) = none_ev              tof12_event(j,i) = none_ev
415              IF ((tof12(j,i,itdc).LT.2000).AND.(tof12(j,i,itdc).GT.100))              IF ((tof12(j,i,itdc).LT.2000).AND.(tof12(j,i,itdc).GT.100))
416       +           tof12_event(j,i) = tof12_event(j,i) + tdc_ev       +           tof12_event(j,i) = tof12_event(j,i) + tdc_ev
             IF ((tof12(j,i,iadc).GT.secure).AND.  
      +           (tof12(j,i,iadc).LT.4095))  
      +           tof12_event(j,i) = tof12_event(j,i) + adc_ev  
417           ENDDO           ENDDO
418        ENDDO        ENDDO
419    
# Line 265  C     middle tof  S21 Line 444  C     middle tof  S21
444              tof21_event(j,i) = none_ev              tof21_event(j,i) = none_ev
445              IF ((tof21(j,i,itdc).LT.2000).AND.(tof21(j,i,itdc).GT.100))              IF ((tof21(j,i,itdc).LT.2000).AND.(tof21(j,i,itdc).GT.100))
446       +           tof21_event(j,i) = tof21_event(j,i) + tdc_ev       +           tof21_event(j,i) = tof21_event(j,i) + tdc_ev
             IF ((tof21(j,i,iadc).GT.secure).AND.  
      +           (tof21(j,i,iadc).LT.4095))  
      +           tof21_event(j,i) = tof21_event(j,i) + adc_ev  
447           ENDDO           ENDDO
448        ENDDO        ENDDO
449    
# Line 297  C     middle tof  S22 Line 473  C     middle tof  S22
473              tof22_event(j,i) = none_ev              tof22_event(j,i) = none_ev
474              IF ((tof22(j,i,itdc).LT.2000).AND.(tof22(j,i,itdc).GT.100))              IF ((tof22(j,i,itdc).LT.2000).AND.(tof22(j,i,itdc).GT.100))
475       +           tof22_event(j,i) = tof22_event(j,i) + tdc_ev       +           tof22_event(j,i) = tof22_event(j,i) + tdc_ev
             IF ((tof22(j,i,iadc).GT.secure).AND.  
      +           (tof22(j,i,iadc).LT.4095))  
      +           tof22_event(j,i) = tof22_event(j,i) + adc_ev  
476           ENDDO           ENDDO
477        ENDDO        ENDDO
478    
# Line 330  C     bottom tof  S31 Line 503  C     bottom tof  S31
503              tof31_event(j,i) = none_ev              tof31_event(j,i) = none_ev
504              IF ((tof31(j,i,itdc).LT.2000).AND.(tof31(j,i,itdc).GT.100))              IF ((tof31(j,i,itdc).LT.2000).AND.(tof31(j,i,itdc).GT.100))
505       +           tof31_event(j,i) = tof31_event(j,i) + tdc_ev       +           tof31_event(j,i) = tof31_event(j,i) + tdc_ev
             IF ((tof31(j,i,iadc).GT.secure).AND.  
      +           (tof31(j,i,iadc).LT.4095))  
      +           tof31_event(j,i) = tof31_event(j,i) + adc_ev  
506           ENDDO           ENDDO
507        ENDDO        ENDDO
508    
# Line 362  C     bottom tof  S32 Line 532  C     bottom tof  S32
532              tof32_event(j,i) = none_ev              tof32_event(j,i) = none_ev
533              IF ((tof32(j,i,itdc).LT.2000).AND.(tof32(j,i,itdc).GT.100))              IF ((tof32(j,i,itdc).LT.2000).AND.(tof32(j,i,itdc).GT.100))
534       +           tof32_event(j,i) = tof32_event(j,i) + tdc_ev       +           tof32_event(j,i) = tof32_event(j,i) + tdc_ev
             IF ((tof32(j,i,iadc).GT.secure).AND.  
      +           (tof32(j,i,iadc).LT.4095))  
      +           tof32_event(j,i) = tof32_event(j,i) + adc_ev  
535           ENDDO           ENDDO
536        ENDDO        ENDDO
537    
# Line 386  c     check if an other paddle has also Line 553  c     check if an other paddle has also
553                 ENDIF                 ENDIF
554              ENDIF              ENDIF
555           ENDIF           ENDIF
556        ENDDO         ENDDO
557    
558        do i=1,6         do i=1,6
559           tof_i_flag(i)=0           tof_i_flag(i)=0
560           tof_j_flag(i)=0           tof_j_flag(i)=0
561        enddo         enddo
562    
563        tof_i_flag(1)=tof11_i         tof_i_flag(1)=tof11_i
564        tof_i_flag(2)=tof12_i         tof_i_flag(2)=tof12_i
565        tof_i_flag(3)=tof21_i         tof_i_flag(3)=tof21_i
566        tof_i_flag(4)=tof22_i         tof_i_flag(4)=tof22_i
567        tof_i_flag(5)=tof31_i         tof_i_flag(5)=tof31_i
568        tof_i_flag(6)=tof32_i         tof_i_flag(6)=tof32_i
569    
570        tof_j_flag(1)=tof11_j         tof_j_flag(1)=tof11_j
571        tof_j_flag(2)=tof12_j         tof_j_flag(2)=tof12_j
572        tof_j_flag(3)=tof21_j         tof_j_flag(3)=tof21_j
573        tof_j_flag(4)=tof22_j         tof_j_flag(4)=tof22_j
574        tof_j_flag(5)=tof31_j         tof_j_flag(5)=tof31_j
575        tof_j_flag(6)=tof32_j         tof_j_flag(6)=tof32_j
576    
577           hitvec(1)=tof11_i
578           hitvec(2)=tof12_i
579           hitvec(3)=tof21_i
580           hitvec(4)=tof22_i
581           hitvec(5)=tof31_i
582           hitvec(6)=tof32_i
583    
584    
 C--------------------------------------------------------------------  
 C--------------------Time walk correction  -------------------------  
 C--------------------------------------------------------------------  
   
       DO i=1,8  
        xhelp_a = tof11(left,i,iadc)  
        xhelp_t = tof11(left,i,itdc)  
        if(xhelp_a>0) xhelp = tw11(left,i)/sqrt(xhelp_a)  
        tof11(left,i,itdc) = xhelp_t  + xhelp  
        tdc_c(ch11a(i),hb11a(i))=tof11(left,i,itdc)  
        xhelp_a = tof11(right,i,iadc)  
        xhelp_t = tof11(right,i,itdc)  
        if(xhelp_a>0) xhelp = tw11(right,i)/sqrt(xhelp_a)  
        tof11(right,i,itdc) = xhelp_t  + xhelp  
        tdc_c(ch11b(i),hb11b(i))=tof11(right,i,itdc)  
       ENDDO  
   
       DO i=1,6  
        xhelp_a = tof12(left,i,iadc)  
        xhelp_t = tof12(left,i,itdc)  
        if(xhelp_a>0) xhelp = tw12(left,i)/sqrt(xhelp_a)  
        tof12(left,i,itdc) = xhelp_t  + xhelp  
        tdc_c(ch12a(i),hb12a(i))=tof12(left,i,itdc)  
        xhelp_a = tof12(right,i,iadc)  
        xhelp_t = tof12(right,i,itdc)  
        if(xhelp_a>0) xhelp = tw12(right,i)/sqrt(xhelp_a)  
        tof12(right,i,itdc) = xhelp_t  + xhelp  
        tdc_c(ch12b(i),hb12b(i))=tof12(right,i,itdc)  
       ENDDO  
 C----  
       DO i=1,2  
        xhelp_a = tof21(left,i,iadc)  
        xhelp_t = tof21(left,i,itdc)  
        if(xhelp_a>0) xhelp = tw21(left,i)/sqrt(xhelp_a)  
        tof21(left,i,itdc) = xhelp_t  + xhelp  
        tdc_c(ch21a(i),hb21a(i))=tof21(left,i,itdc)  
        xhelp_a = tof21(right,i,iadc)  
        xhelp_t = tof21(right,i,itdc)  
        if(xhelp_a>0) xhelp = tw21(right,i)/sqrt(xhelp_a)  
        tof21(right,i,itdc) = xhelp_t  + xhelp  
        tdc_c(ch21b(i),hb21b(i))=tof21(right,i,itdc)  
       ENDDO  
   
       DO i=1,2  
        xhelp_a = tof22(left,i,iadc)  
        xhelp_t = tof22(left,i,itdc)  
        if(xhelp_a>0) xhelp = tw22(left,i)/sqrt(xhelp_a)  
        tof22(left,i,itdc) = xhelp_t  + xhelp  
        tdc_c(ch22a(i),hb22a(i))=tof22(left,i,itdc)  
        xhelp_a = tof22(right,i,iadc)  
        xhelp_t = tof22(right,i,itdc)  
        if(xhelp_a>0) xhelp = tw22(right,i)/sqrt(xhelp_a)  
        tof22(right,i,itdc) = xhelp_t  + xhelp  
        tdc_c(ch22b(i),hb22b(i))=tof22(right,i,itdc)  
       ENDDO  
 C----  
   
       DO i=1,3  
        xhelp_a = tof31(left,i,iadc)  
        xhelp_t = tof31(left,i,itdc)  
        if(xhelp_a>0) xhelp = tw31(left,i)/sqrt(xhelp_a)  
        tof31(left,i,itdc) = xhelp_t  + xhelp  
        tdc_c(ch31a(i),hb31a(i))=tof31(left,i,itdc)  
        xhelp_a = tof31(right,i,iadc)  
        xhelp_t = tof31(right,i,itdc)  
        if(xhelp_a>0) xhelp = tw31(right,i)/sqrt(xhelp_a)  
        tof31(right,i,itdc) = xhelp_t  + xhelp  
        tdc_c(ch31b(i),hb31b(i))=tof31(right,i,itdc)  
       ENDDO  
   
       DO i=1,3  
        xhelp_a = tof32(left,i,iadc)  
        xhelp_t = tof32(left,i,itdc)  
        if(xhelp_a>0) xhelp = tw32(left,i)/sqrt(xhelp_a)  
        tof32(left,i,itdc) = xhelp_t  + xhelp  
        tdc_c(ch32a(i),hb32a(i))=tof32(left,i,itdc)  
        xhelp_a = tof32(right,i,iadc)  
        xhelp_t = tof32(right,i,itdc)  
        if(xhelp_a>0) xhelp = tw32(right,i)/sqrt(xhelp_a)  
        tof32(right,i,itdc) = xhelp_t  + xhelp  
        tdc_c(ch32b(i),hb32b(i))=tof32(right,i,itdc)  
       ENDDO  
 C----  
   
585  C------------------------------------------------------------------  C------------------------------------------------------------------
586  C---  calculate track position in paddle using timing difference  C--  calculate track position in paddle using timing difference
587    C--  this calculation is preliminary and uses some standard
588    C--  calibration values, but we need to find a rough position to
589    C--  be able to calculate artificial ADC values (needed for the
590    C--  timewalk...
591  C------------------------------------------------------------------  C------------------------------------------------------------------
592    
593        do i=1,3         do i=1,3
594           xtofpos(i)=100.           xtofpos(i)=100.
595           ytofpos(i)=100.           ytofpos(i)=100.
596        enddo         enddo
597    
598  C-----------------------------S1 --------------------------------  C-----------------------------S1 --------------------------------
599    
600        IF (tof11_i.GT.none_find) THEN        IF (tof11_i.GT.none_find) THEN
601           ytofpos(1)  = ((tof11(1,tof11_i,itdc)-tof11(2,tof11_i,itdc))/2.           ytofpos(1)  = ((tof11(1,tof11_i,itdc)-tof11(2,tof11_i,itdc))/2.
602       +        -y_coor_lin11(tof11_i,offset))/y_coor_lin11(tof11_i,slope)       +   -y_coor_lin11c(tof11_i,offset))/y_coor_lin11c(tof11_i,slope)
603        endif        endif
604    
605        IF (tof12_i.GT.none_find) THEN        IF (tof12_i.GT.none_find) THEN
606           xtofpos(1)  = ((tof12(1,tof12_i,itdc)-tof12(2,tof12_i,itdc))/2.           xtofpos(1)  = ((tof12(1,tof12_i,itdc)-tof12(2,tof12_i,itdc))/2.
607       +        -x_coor_lin12(tof12_i,offset))/x_coor_lin12(tof12_i,slope)       +   -x_coor_lin12c(tof12_i,offset))/x_coor_lin12c(tof12_i,slope)
608        endif        endif
609    
610    
# Line 517  C-----------------------------S2 ------- Line 612  C-----------------------------S2 -------
612    
613        IF (tof21_i.GT.none_find) THEN        IF (tof21_i.GT.none_find) THEN
614           xtofpos(2) = ((tof21(1,tof21_i,itdc)-tof21(2,tof21_i,itdc))/2.           xtofpos(2) = ((tof21(1,tof21_i,itdc)-tof21(2,tof21_i,itdc))/2.
615       +        -x_coor_lin21(tof21_i,offset))/x_coor_lin21(tof21_i,slope)       +    -x_coor_lin21c(tof21_i,offset))/x_coor_lin21c(tof21_i,slope)
616        endif        endif
617    
618        IF (tof22_i.GT.none_find) THEN        IF (tof22_i.GT.none_find) THEN
619           ytofpos(2) = ((tof22(1,tof22_i,itdc)-tof22(2,tof22_i,itdc))/2.           ytofpos(2) = ((tof22(1,tof22_i,itdc)-tof22(2,tof22_i,itdc))/2.
620       +        -y_coor_lin22(tof22_i,offset))/y_coor_lin22(tof22_i,slope)       +    -y_coor_lin22c(tof22_i,offset))/y_coor_lin22c(tof22_i,slope)
621        endif        endif
622                
623    
# Line 530  C-----------------------------S3 ------- Line 625  C-----------------------------S3 -------
625    
626        IF (tof31_i.GT.none_find) THEN        IF (tof31_i.GT.none_find) THEN
627           ytofpos(3)  = ((tof31(1,tof31_i,itdc)-tof31(2,tof31_i,itdc))/2.           ytofpos(3)  = ((tof31(1,tof31_i,itdc)-tof31(2,tof31_i,itdc))/2.
628       +        -y_coor_lin31(tof31_i,offset))/y_coor_lin31(tof31_i,slope)       +    -y_coor_lin31c(tof31_i,offset))/y_coor_lin31c(tof31_i,slope)
629        endif        endif
630    
631        IF (tof32_i.GT.none_find) THEN        IF (tof32_i.GT.none_find) THEN
632           xtofpos(3)  = ((tof32(1,tof32_i,itdc)-tof32(2,tof32_i,itdc))/2.           xtofpos(3)  = ((tof32(1,tof32_i,itdc)-tof32(2,tof32_i,itdc))/2.
633       +        -x_coor_lin32(tof32_i,offset))/x_coor_lin32(tof32_i,slope)       +    -x_coor_lin32c(tof32_i,offset))/x_coor_lin32c(tof32_i,slope)
634        endif        endif
635    
636    
       do i=1,3  
          if (abs(xtofpos(i)).gt.100.) then  
             xtofpos(i)=101.  
          endif  
          if (abs(ytofpos(i)).gt.100.) then  
             ytofpos(i)=101.  
          endif  
       enddo  
   
637  C----------------------------------------------------------------------  C----------------------------------------------------------------------
638  C--------------------Corrections on ADC-data -------------------------  C---------------------  zenith angle theta  ---------------------------
 C---------------------zenith angle theta  ---------------------------  
639  C----------------------------------------------------------------------  C----------------------------------------------------------------------
640    
641        dx=0.         xhelp1=0.
642        dy=0.         if (tof11_i.GT.none_find) xhelp1=tof11_x(tof11_i)
643        dr=0.         if (xtofpos(1).lt.100)  xhelp1=xtofpos(1)
644        theta13 = 0.  
645           yhelp1=0.
646           IF ((tof12_i.GT.none_find).AND.(tof32_i.GT.none_find))         if (tof12_i.GT.none_find) yhelp1=tof12_y(tof12_i)
647       &        dx  = xtofpos(1) - xtofpos(3)         if (ytofpos(1).lt.100)  yhelp1=ytofpos(1)
648           IF ((tof11_i.GT.none_find).AND.(tof31_i.GT.none_find))  
649       &        dy  = ytofpos(1) - ytofpos(3)  
650           dr = sqrt(dx*dx+dy*dy)         yhelp2=0.
651           theta13 = atan(dr/tofarm13)         if (tof32_i.GT.none_find) yhelp2=tof32_y(tof32_i)
652           if (ytofpos(3).lt.100)  yhelp2=ytofpos(3)
653        dx=0.  
654        dy=0.         xhelp2=0.
655        dr=0.         if (tof31_i.GT.none_find) xhelp2=tof31_x(tof31_i)
656        theta12 = 0.         if (xtofpos(3).lt.100)  xhelp2=xtofpos(3)
657    
658           IF ((tof12_i.GT.none_find).AND.(tof21_i.GT.none_find))  
659       &        dx  = xtofpos(1) - xtofpos(2)         dx=0.
660           IF ((tof11_i.GT.none_find).AND.(tof22_i.GT.none_find))         dy=0.
661       &        dy  = ytofpos(1) - ytofpos(2)         dr=0.
662           dr = sqrt(dx*dx+dy*dy)         theta13 = 0.
663           theta12 = atan(dr/tofarm12)  
664                   dx  = xhelp1 - xhelp2
665        dx=0.         dy  = yhelp1 - yhelp2
666        dy=0.         dr = sqrt(dx*dx+dy*dy)
667        dr=0.         theta13 = atan(dr/tofarm13)
668        theta23 = 0.  
669    
          IF ((tof21_i.GT.none_find).AND.(tof32_i.GT.none_find))  
      &        dx  = xtofpos(2) - xtofpos(3)  
          IF ((tof22_i.GT.none_find).AND.(tof31_i.GT.none_find))  
      &        dy  = ytofpos(2) - ytofpos(3)  
          dr = sqrt(dx*dx+dy*dy)  
          theta23 = atan(dr/tofarm23)  
           
           
670  C----------------------------------------------------------------------  C----------------------------------------------------------------------
671  C------------------angle and ADC(x) correction  C--- check charge:
672    C--- if Z=2 we should use the attenuation curve for helium to
673    C--- fill the artificail ADC values and NOT divide by "hepratio"
674    C--- if Z>2 we should do a correction to
675    C--- the k1 constants in the beta calculation
676  C----------------------------------------------------------------------  C----------------------------------------------------------------------
677  C-----------------------------S1 --------------------------------  
678             iz = int(check_charge(theta13,hitvec))
679    c         write(*,*) 'charge in tofl2com',iz
680    
681    C--------------------------------------------------------------------
682    C---- if TDCleft.and.TDCright and NO ADC insert artificial ADC
683    C---- values
684    C--------------------------------------------------------------------
685  c middle y (or x) position of the upper and middle ToF-Paddle  c middle y (or x) position of the upper and middle ToF-Paddle
686  c       DATA tof11_x/ -17.85,-12.75,-7.65,-2.55,2.55,7.65,12.75,17.85/  c       DATA tof11_x/ -17.85,-12.75,-7.65,-2.55,2.55,7.65,12.75,17.85/
687  c       DATA tof12_y/ -13.75,-8.25,-2.75,2.75,8.25,13.75/  c       DATA tof12_y/ -13.75,-8.25,-2.75,2.75,8.25,13.75/
688  c       DATA tof21_y/ -3.75,3.75/  c       DATA tof21_y/  3.75,-3.75/     ! paddles in different order
689  c       DATA tof22_x/ -4.5,4.5/  c       DATA tof22_x/ -4.5,4.5/
690  c       DATA tof31_x/ -6.0,0.,6.0/  c       DATA tof31_x/ -6.0,0.,6.0/
691  c       DATA tof32_y/ -5.0,0.0,5.0/  c       DATA tof32_y/ -5.0,0.0,5.0/
692    
       yhelp=0.  
       if (tof12_i.GT.none_find) yhelp=tof12_y(tof12_i)  
       if (ytofpos(1).lt.100)  yhelp=ytofpos(1)  
693    
694        IF (tof11_i.GT.none_find.AND.abs(yhelp).lt.100) THEN  C----------------------------  S1 -------------------------------------
695    
696    c       yhelp=0.
697           yhelp=100.  ! WM
698           if (tof12_i.GT.none_find) yhelp=tof12_y(tof12_i)
699           if (ytofpos(1).lt.100)  yhelp=ytofpos(1)
700    
701           IF (tof11_i.GT.none_find.AND.abs(yhelp).lt.100) THEN
702           i = tof11_i           i = tof11_i
703           xdummy=tof11(left,i,iadc)           if (adc(ch11a(i),hb11a(i)).eq.4095) then
704           tof11(left,i,iadc) = tof11(left,i,iadc)*cos(theta13)              xkorr = atten(left,11,i,yhelp)
705           if (tof11(left,i,iadc).lt.1000) then              if (iz.le.1) xkorr=xkorr/hepratio
706              xkorr=adcx11(left,i,1)*exp(-yhelp/adcx11(left,i,2))              tof11(left,i,iadc)=xkorr/cos(theta13)
707              xkorr0=adcx11(left,i,1)              adcflagtof(ch11a(i),hb11a(i)) = 1
             adctof_c(ch11a(i),hb11a(i))=tof11(left,i,iadc)/xkorr  
708           endif           endif
709             if (adc(ch11b(i),hb11b(i)).eq.4095) then
710           tof11(right,i,iadc) = tof11(right,i,iadc)*cos(theta13)              xkorr = atten(right,11,i,yhelp)
711           if (tof11(right,i,iadc).lt.1000) then              if (iz.le.1) xkorr=xkorr/hepratio
712              xkorr=adcx11(right,i,1)*exp(yhelp/adcx11(right,i,2))              tof11(right,i,iadc)=xkorr/cos(theta13)
713              xkorr0=adcx11(right,i,1)              adcflagtof(ch11b(i),hb11b(i)) = 1
             adctof_c(ch11b(i),hb11b(i))=tof11(right,i,iadc)/xkorr  
714           endif           endif
715        ENDIF         ENDIF
   
       xhelp=0.  
       if (tof11_i.GT.none_find) xhelp=tof11_x(tof11_i)  
       if (xtofpos(1).lt.100)  xhelp=xtofpos(1)  
716    
717        IF (tof12_i.GT.none_find.AND.abs(xhelp).lt.100) THEN  c       xhelp=0.
718           xhelp=100.  ! WM
719           if (tof11_i.GT.none_find) xhelp=tof11_x(tof11_i)
720           if (xtofpos(1).lt.100)  xhelp=xtofpos(1)
721    
722           IF (tof12_i.GT.none_find.AND.abs(xhelp).lt.100) THEN
723           i = tof12_i           i = tof12_i
724           tof12(left,i,iadc) = tof12(left,i,iadc)*cos(theta13)           if (adc(ch12a(i),hb12a(i)).eq.4095) then
725           if (tof12(left,i,iadc).lt.1000) then              xkorr = atten(left,12,i,xhelp)
726              xkorr=adcx12(left,i,1)*exp(-xhelp/adcx12(left,i,2))              if (iz.le.1) xkorr=xkorr/hepratio
727              xkorr0=adcx12(left,i,1)              tof12(left,i,iadc) = xkorr/cos(theta13)
728              adctof_c(ch12a(i),hb12a(i))=tof12(left,i,iadc)/xkorr              adcflagtof(ch12a(i),hb12a(i)) = 1
729           endif           endif
730             if (adc(ch12b(i),hb12b(i)).eq.4095) then
731           tof12(right,i,iadc) = tof12(right,i,iadc)*cos(theta13)              xkorr = atten(right,12,i,xhelp)
732           if (tof12(right,i,iadc).lt.1000) then              if (iz.le.1) xkorr=xkorr/hepratio
733              xkorr=adcx12(right,i,1)*exp(xhelp/adcx12(right,i,2))              tof12(right,i,iadc) = xkorr/cos(theta13)
734              xkorr0=adcx12(right,i,1)              adcflagtof(ch12b(i),hb12b(i)) = 1
             adctof_c(ch12b(i),hb12b(i))=tof12(right,i,iadc)/xkorr  
735           endif           endif
736        ENDIF         ENDIF
737    
738  C-----------------------------S2 --------------------------------  C-----------------------------S2 --------------------------------
739    
740        xhelp=0.  c       xhelp=0.
741        if (tof22_i.GT.none_find) xhelp=tof22_x(tof22_i)         xhelp=100.   ! WM
742        if (xtofpos(2).lt.100)  xhelp=xtofpos(2)         if (tof22_i.GT.none_find) xhelp=tof22_x(tof22_i)
743           if (xtofpos(2).lt.100)  xhelp=xtofpos(2)
       IF (tof21_i.GT.none_find.AND.abs(xhelp).lt.100) THEN  
744    
745           IF (tof21_i.GT.none_find.AND.abs(xhelp).lt.100) THEN
746           i = tof21_i           i = tof21_i
747           tof21(left,i,iadc) = tof21(left,i,iadc)*cos(theta13)           if (adc(ch21a(i),hb21a(i)).eq.4095) then
748           if (tof21(left,i,iadc).lt.1000) then              xkorr = atten(left,21,i,xhelp)
749              xkorr=adcx21(left,i,1)*exp(-xhelp/adcx21(left,i,2))              if (iz.le.1) xkorr=xkorr/hepratio
750              xkorr0=adcx21(left,i,1)              tof21(left,i,iadc) = xkorr/cos(theta13)
751              adctof_c(ch21a(i),hb21a(i))=tof21(left,i,iadc)/xkorr              adcflagtof(ch21a(i),hb21a(i)) = 1
752           endif           endif
753             if (adc(ch21b(i),hb21b(i)).eq.4095) then
754           tof21(right,i,iadc) = tof21(right,i,iadc)*cos(theta13)              xkorr = atten(right,21,i,xhelp)
755           if (tof21(right,i,iadc).lt.1000) then              if (iz.le.1) xkorr=xkorr/hepratio
756              xkorr=adcx21(right,i,1)*exp(xhelp/adcx21(right,i,2))              tof21(right,i,iadc) = xkorr/cos(theta13)
757              xkorr0=adcx21(right,i,1)              adcflagtof(ch21b(i),hb21b(i)) = 1
             adctof_c(ch21b(i),hb21b(i))=tof21(right,i,iadc)/xkorr  
758           endif           endif
759        ENDIF         ENDIF
   
760    
       yhelp=0.  
       if (tof21_i.GT.none_find) yhelp=tof21_y(tof21_i)  
       if (ytofpos(2).lt.100)  yhelp=ytofpos(2)  
761    
762        IF (tof22_i.GT.none_find.AND.abs(yhelp).lt.100) THEN  c       yhelp=0.
763           yhelp=100.   ! WM
764           if (tof21_i.GT.none_find) yhelp=tof21_y(tof21_i)
765           if (ytofpos(2).lt.100)  yhelp=ytofpos(2)
766    
767           IF (tof22_i.GT.none_find.AND.abs(yhelp).lt.100) THEN
768           i = tof22_i           i = tof22_i
769           tof22(left,i,iadc) = tof22(left,i,iadc)*cos(theta13)           if (adc(ch22a(i),hb22a(i)).eq.4095) then
770           if (tof22(left,i,iadc).lt.1000) then              xkorr = atten(left,22,i,yhelp)
771              xkorr=adcx22(left,i,1)*exp(-yhelp/adcx22(left,i,2))              if (iz.le.1) xkorr=xkorr/hepratio
772              xkorr0=adcx22(left,i,1)              tof22(left,i,iadc) = xkorr/cos(theta13)
773              adctof_c(ch22a(i),hb22a(i))=tof22(left,i,iadc)/xkorr              adcflagtof(ch22a(i),hb22a(i)) = 1
774           endif           endif
775             if (adc(ch22b(i),hb22b(i)).eq.4095) then
776           tof22(right,i,iadc) = tof22(right,i,iadc)*cos(theta13)              xkorr = atten(right,22,i,yhelp)
777           if (tof22(right,i,iadc).lt.1000) then              if (iz.le.1) xkorr=xkorr/hepratio
778              xkorr=adcx22(right,i,1)*exp(yhelp/adcx22(right,i,2))              tof22(right,i,iadc) = xkorr/cos(theta13)
779              xkorr0=adcx22(right,i,1)              adcflagtof(ch22b(i),hb22b(i)) = 1
             adctof_c(ch22b(i),hb22b(i))=tof22(right,i,iadc)/xkorr  
780           endif           endif
781        ENDIF         ENDIF
782    
783  C-----------------------------S3 --------------------------------  C-----------------------------S3 --------------------------------
784    
785        yhelp=0.  c       yhelp=0.
786        if (tof32_i.GT.none_find) yhelp=tof32_y(tof32_i)         yhelp=100.  ! WM
787        if (ytofpos(3).lt.100)  yhelp=ytofpos(3)         if (tof32_i.GT.none_find) yhelp=tof32_y(tof32_i)
788           if (ytofpos(3).lt.100)  yhelp=ytofpos(3)
       IF (tof31_i.GT.none_find.AND.abs(yhelp).lt.100) THEN  
789    
790           IF (tof31_i.GT.none_find.AND.abs(yhelp).lt.100) THEN
791           i = tof31_i           i = tof31_i
792           tof31(left,i,iadc) = tof31(left,i,iadc)*cos(theta13)           if (adc(ch31a(i),hb31a(i)).eq.4095) then
793           if (tof31(left,i,iadc).lt.1000) then              xkorr = atten(left,31,i,yhelp)
794              xkorr=adcx31(left,i,1)*exp(-yhelp/adcx31(left,i,2))              if (iz.le.1) xkorr=xkorr/hepratio
795              xkorr0=adcx31(left,i,1)              tof31(left,i,iadc) = xkorr/cos(theta13)
796              adctof_c(ch31a(i),hb31a(i))=tof31(left,i,iadc)/xkorr              adcflagtof(ch31a(i),hb31a(i)) = 1
797           endif           endif
798             if (adc(ch31b(i),hb31b(i)).eq.4095) then
799           tof31(right,i,iadc) = tof31(right,i,iadc)*cos(theta13)              xkorr = atten(right,31,i,yhelp)
800           if (tof31(right,i,iadc).lt.1000) then              if (iz.le.1) xkorr=xkorr/hepratio
801              xkorr=adcx31(right,i,1)*exp(yhelp/adcx31(right,i,2))              tof31(right,i,iadc) = xkorr/cos(theta13)
802              xkorr0=adcx31(right,i,1)              adcflagtof(ch31b(i),hb31b(i)) = 1
             adctof_c(ch31b(i),hb31b(i))=tof31(right,i,iadc)/xkorr  
803           endif           endif
804        ENDIF         ENDIF
805    
806        xhelp=0.  c       xhelp=0.
807        if (tof31_i.GT.none_find) xhelp=tof31_x(tof31_i)         xhelp=100.   ! WM
808        if (xtofpos(3).lt.100)  xhelp=xtofpos(3)         if (tof31_i.GT.none_find) xhelp=tof31_x(tof31_i)
809           if (xtofpos(3).lt.100)  xhelp=xtofpos(3)
       IF (tof32_i.GT.none_find.AND.abs(xhelp).lt.100) THEN  
810    
811           IF (tof32_i.GT.none_find.AND.abs(xhelp).lt.100) THEN
812           i = tof32_i           i = tof32_i
813           tof32(left,i,iadc) = tof32(left,i,iadc)*cos(theta13)           if (adc(ch32a(i),hb32a(i)).eq.4095) then
814           if (tof32(left,i,iadc).lt.1000) then              xkorr = atten(left,32,i,xhelp)
815              xkorr=adcx32(left,i,1)*exp(-xhelp/adcx32(left,i,2))              if (iz.le.1) xkorr=xkorr/hepratio
816              xkorr0=adcx32(left,i,1)              tof32(left,i,iadc) = xkorr/cos(theta13)
817              adctof_c(ch32a(i),hb32a(i))=tof32(left,i,iadc)/xkorr              adcflagtof(ch32a(i),hb32a(i)) = 1
818           endif           endif
819             if (adc(ch32b(i),hb32b(i)).eq.4095) then
820           tof32(right,i,iadc) = tof32(right,i,iadc)*cos(theta13)              xkorr = atten(right,32,i,xhelp)
821           if (tof32(right,i,iadc).lt.1000) then              if (iz.le.1) xkorr=xkorr/hepratio
822              xkorr=adcx32(right,i,1)*exp(xhelp/adcx32(right,i,2))              tof32(right,i,iadc) = xkorr/cos(theta13)
823              xkorr0=adcx32(right,i,1)              adcflagtof(ch32b(i),hb32b(i)) = 1
             adctof_c(ch32b(i),hb32b(i))=tof32(right,i,iadc)/xkorr  
824           endif           endif
825        ENDIF         ENDIF
826    
827    
828    C-------------------------------------------------------------------
829    C--------------------Time walk correction  -------------------------
830    C-------------------------------------------------------------------
831    C-------------------------------------------------------------------
832    C Now there is for each hitted paddle a TDC and ADC value, if the
833    C TDC was < 4095.
834    C There might be also TDC-ADC pairs in paddles not hitted
835    
836    C-------------------------------------------------------------------
837    C If we have multiple paddles hit, so that no artificial ADC value
838    C is created, we set the raw TDC value as "tdc_c"
839    C-------------------------------------------------------------------
840    c
841    c       do i=1,4
842    c         do j=1,12
843    c            tdc_c(i,j) = tdc(i,j)
844    c         enddo
845    c       enddo
846    c
847    C----  Let's correct the raw TDC value with the time walk  ---------
848    
849          DO i=1,8
850             if ((tdc(ch11a(i),hb11a(i)).lt.4095).and.
851         &             (tof11(left,i,iadc).lt.3786)) THEN
852             xhelp = tw11(left,i)/(tof11(left,i,iadc)**0.5)
853             tof11(left,i,itdc) = tof11(left,i,itdc) + xhelp
854             tdc_c(ch11a(i),hb11a(i))=tof11(left,i,itdc)
855                                                  ENDIF
856    
857             if ((tdc(ch11b(i),hb11b(i)).lt.4095).and.
858         &             (tof11(right,i,iadc).lt.3786)) THEN
859             xhelp = tw11(right,i)/(tof11(right,i,iadc)**0.5)
860             tof11(right,i,itdc) = tof11(right,i,itdc) + xhelp
861             tdc_c(ch11b(i),hb11b(i))=tof11(right,i,itdc)
862                                                 ENDIF
863          ENDDO
864    
865    
866          DO i=1,6
867             if ((tdc(ch12a(i),hb12a(i)).lt.4095).and.
868         &             (tof12(left,i,iadc).lt.3786)) THEN
869             xhelp = tw12(left,i)/(tof12(left,i,iadc)**0.5)
870             tof12(left,i,itdc) = tof12(left,i,itdc) + xhelp
871             tdc_c(ch12a(i),hb12a(i))=tof12(left,i,itdc)
872                                                  ENDIF
873    
874             if ((tdc(ch12b(i),hb12b(i)).lt.4095).and.
875         &             (tof12(right,i,iadc).lt.3786)) THEN
876             xhelp = tw12(right,i)/(tof12(right,i,iadc)**0.5)
877             tof12(right,i,itdc) = tof12(right,i,itdc) + xhelp
878             tdc_c(ch12b(i),hb12b(i))=tof12(right,i,itdc)
879                                                 ENDIF
880          ENDDO
881    
882    C----
883          DO I=1,2
884             if ((tdc(ch21a(i),hb21a(i)).lt.4095).and.
885         &             (tof21(left,i,iadc).lt.3786)) THEN
886             xhelp = tw21(left,i)/(tof21(left,i,iadc)**0.5)
887             tof21(left,i,itdc) = tof21(left,i,itdc) + xhelp
888             tdc_c(ch21a(i),hb21a(i))=tof21(left,i,itdc)
889                                                  ENDIF
890    
891             if ((tdc(ch21b(i),hb21b(i)).lt.4095).and.
892         &             (tof21(right,i,iadc).lt.3786)) THEN
893             xhelp = tw21(right,i)/(tof21(right,i,iadc)**0.5)
894             tof21(right,i,itdc) = tof21(right,i,itdc) + xhelp
895             tdc_c(ch21b(i),hb21b(i))=tof21(right,i,itdc)
896                                                 ENDIF
897          ENDDO
898    
899          DO I=1,2
900             if ((tdc(ch22a(i),hb22a(i)).lt.4095).and.
901         &             (tof22(left,i,iadc).lt.3786)) THEN
902             xhelp = tw22(left,i)/(tof22(left,i,iadc)**0.5)
903             tof22(left,i,itdc) = tof22(left,i,itdc) + xhelp
904             tdc_c(ch22a(i),hb22a(i))=tof22(left,i,itdc)
905                                                  ENDIF
906    
907             if ((tdc(ch22b(i),hb22b(i)).lt.4095).and.
908         &             (tof22(right,i,iadc).lt.3786)) THEN
909             xhelp = tw22(right,i)/(tof22(right,i,iadc)**0.5)
910             tof22(right,i,itdc) = tof22(right,i,itdc) + xhelp
911             tdc_c(ch22b(i),hb22b(i))=tof22(right,i,itdc)
912                                                 ENDIF
913          ENDDO
914    
915    C----
916          DO I=1,3
917             if ((tdc(ch31a(i),hb31a(i)).lt.4095).and.
918         &             (tof31(left,i,iadc).lt.3786)) THEN
919             xhelp = tw31(left,i)/(tof31(left,i,iadc)**0.5)
920             tof31(left,i,itdc) = tof31(left,i,itdc) + xhelp
921             tdc_c(ch31a(i),hb31a(i))=tof31(left,i,itdc)
922                                                  ENDIF
923    
924             if ((tdc(ch31b(i),hb31b(i)).lt.4095).and.
925         &             (tof31(right,i,iadc).lt.3786)) THEN
926             xhelp = tw31(right,i)/(tof31(right,i,iadc)**0.5)
927             tof31(right,i,itdc) = tof31(right,i,itdc) + xhelp
928             tdc_c(ch31b(i),hb31b(i))=tof31(right,i,itdc)
929                                                 ENDIF
930          ENDDO
931    
932          DO I=1,3
933             if ((tdc(ch32a(i),hb32a(i)).lt.4095).and.
934         &             (tof32(left,i,iadc).lt.3786)) THEN
935             xhelp = tw32(left,i)/(tof32(left,i,iadc)**0.5)
936             tof32(left,i,itdc) = tof32(left,i,itdc) + xhelp
937             tdc_c(ch32a(i),hb32a(i))=tof32(left,i,itdc)
938                                                  ENDIF
939    
940             if ((tdc(ch32b(i),hb32b(i)).lt.4095).and.
941         &             (tof32(right,i,iadc).lt.3786)) THEN
942             xhelp = tw32(right,i)/(tof32(right,i,iadc)**0.5)
943             tof32(right,i,itdc) = tof32(right,i,itdc) + xhelp
944             tdc_c(ch32b(i),hb32b(i))=tof32(right,i,itdc)
945                                                 ENDIF
946          ENDDO
947          
948    
949    C---------------------------------------------------------------
950    C--- calculate track position in paddle using timing difference
951    C--- now using the time-walk corrected TDC values
952    C---------------------------------------------------------------
953    
954          do i=1,3
955             xtofpos(i)=100.
956             ytofpos(i)=100.
957          enddo
958    
959    C-----------------------------S1 --------------------------------
960    
961          IF (tof11_i.GT.none_find) THEN
962             ytofpos(1)  = ((tof11(1,tof11_i,itdc)-tof11(2,tof11_i,itdc))/2.
963         +        -y_coor_lin11(tof11_i,offset))/y_coor_lin11(tof11_i,slope)
964           i=tof11_i
965          endif
966    
967          IF (tof12_i.GT.none_find) THEN
968             xtofpos(1)  = ((tof12(1,tof12_i,itdc)-tof12(2,tof12_i,itdc))/2.
969         +        -x_coor_lin12(tof12_i,offset))/x_coor_lin12(tof12_i,slope)
970          i=tof12_i
971          endif
972    
973    
974    C-----------------------------S2 --------------------------------
975    
976          IF (tof21_i.GT.none_find) THEN
977             xtofpos(2) = ((tof21(1,tof21_i,itdc)-tof21(2,tof21_i,itdc))/2.
978         +        -x_coor_lin21(tof21_i,offset))/x_coor_lin21(tof21_i,slope)
979          i=tof21_i
980          endif
981    
982          IF (tof22_i.GT.none_find) THEN
983             ytofpos(2) = ((tof22(1,tof22_i,itdc)-tof22(2,tof22_i,itdc))/2.
984         +        -y_coor_lin22(tof22_i,offset))/y_coor_lin22(tof22_i,slope)
985          i=tof22_i
986          endif
987          
988    
989    C-----------------------------S3 --------------------------------
990    
991          IF (tof31_i.GT.none_find) THEN
992             ytofpos(3)  = ((tof31(1,tof31_i,itdc)-tof31(2,tof31_i,itdc))/2.
993         +        -y_coor_lin31(tof31_i,offset))/y_coor_lin31(tof31_i,slope)
994          i=tof31_i
995          endif
996    
997          IF (tof32_i.GT.none_find) THEN
998             xtofpos(3)  = ((tof32(1,tof32_i,itdc)-tof32(2,tof32_i,itdc))/2.
999         +        -x_coor_lin32(tof32_i,offset))/x_coor_lin32(tof32_i,slope)
1000          i=tof32_i
1001          endif
1002    
1003    
1004    c      do i=1,3
1005    c         if (abs(xtofpos(i)).gt.100.) then
1006    c            xtofpos(i)=101.
1007    c         endif
1008    c         if (abs(ytofpos(i)).gt.100.) then
1009    c            ytofpos(i)=101.
1010    c         endif
1011    c      enddo
1012    
1013    
1014    C--  restrict TDC measurements to physical paddle dimensions +/- 10 cm
1015    C--  this cut is now stronger than in the old versions
1016    
1017            if (abs(xtofpos(1)).gt.31.)  xtofpos(1)=101.
1018            if (abs(xtofpos(2)).gt.19.)  xtofpos(2)=101.
1019            if (abs(xtofpos(3)).gt.19.)  xtofpos(3)=101.
1020    
1021            if (abs(ytofpos(1)).gt.26.)  ytofpos(1)=101.
1022            if (abs(ytofpos(2)).gt.18.)  ytofpos(2)=101.
1023            if (abs(ytofpos(3)).gt.18.)  ytofpos(3)=101.
1024    
1025    C----------------------------------------------------------------------
1026    C---------------------  zenith angle theta  ---------------------------
1027    C----------------------------------------------------------------------
1028    C-------------------  improved calculation  ---------------------------
1029    
1030           xhelp1=0.
1031           if (tof11_i.GT.none_find) xhelp1=tof11_x(tof11_i)
1032           if (xtofpos(1).lt.100)  xhelp1=xtofpos(1)
1033    
1034           yhelp1=0.
1035           if (tof12_i.GT.none_find) yhelp1=tof12_y(tof12_i)
1036           if (ytofpos(1).lt.100)  yhelp1=ytofpos(1)
1037    
1038           yhelp2=0.
1039           if (tof32_i.GT.none_find) yhelp2=tof32_y(tof32_i)
1040           if (ytofpos(3).lt.100)  yhelp2=ytofpos(3)
1041    
1042           xhelp2=0.
1043           if (tof31_i.GT.none_find) xhelp2=tof31_x(tof31_i)
1044           if (xtofpos(3).lt.100)  xhelp2=xtofpos(3)
1045    
1046    
1047           dx=0.
1048           dy=0.
1049           dr=0.
1050           theta13 = 0.
1051    
1052           dx  = xhelp1 - xhelp2
1053           dy  = yhelp1 - yhelp2
1054           dr = sqrt(dx*dx+dy*dy)
1055           theta13 = atan(dr/tofarm13)
1056    
1057  C-----------------------------------------------------------------------  
1058    C------------------------------------------------------------------
1059    C------------------------------------------------------------------
1060    C-------angle and ADC(x) correction: moved to new dEdx routine
1061    C------------------------------------------------------------------
1062    C------------------------------------------------------------------
1063    
1064    C--------------------------------------------------------------------
1065  C----------------------calculate Beta  ------------------------------  C----------------------calculate Beta  ------------------------------
1066  C-----------------------------------------------------------------------  C--------------------------------------------------------------------
1067  C-------------------difference of sums  ---------------------------  C-------------------difference of sums  -----------------------------
1068  C  C
1069  C     DS = (t1+t2) - t3+t4)  C     DS = (t1+t2) - t3+t4)
1070  C     DS = c1 + c2/beta*cos(theta)  C     DS = c1 + c2/beta*cos(theta)
# Line 752  C     =>  c2 =  ca.60 for  0.45 m    c2 Line 1073  C     =>  c2 =  ca.60 for  0.45 m    c2
1073  C     since TDC resolution varies slightly c2 has to be calibrated  C     since TDC resolution varies slightly c2 has to be calibrated
1074    
1075  C     S11 - S31  C     S11 - S31
1076        IF (tof11_i.GT.none_find.AND.tof31_i.GT.none_find) THEN  
1077           IF ((tof11_i.GT.none_find).AND.(tof31_i.GT.none_find).AND.
1078         &    (ytofpos(1).NE.101.).AND.(ytofpos(3).NE.101.)) THEN
1079           xhelp1 = tof11(1,tof11_i,itdc)+tof11(2,tof11_i,itdc)           xhelp1 = tof11(1,tof11_i,itdc)+tof11(2,tof11_i,itdc)
1080           xhelp2 = tof31(1,tof31_i,itdc)+tof31(2,tof31_i,itdc)           xhelp2 = tof31(1,tof31_i,itdc)+tof31(2,tof31_i,itdc)
1081           ds = xhelp1-xhelp2           ds = xhelp1-xhelp2
1082           ihelp=(tof11_i-1)*3+tof31_i           ihelp=(tof11_i-1)*3+tof31_i
1083           c1 = k_S11S31(1,ihelp)           c1 = k_S11S31(1,ihelp)
1084             if (iz.gt.2) c1 = c1 + k1corrA1
1085           c2 = k_S11S31(2,ihelp)           c2 = k_S11S31(2,ihelp)
1086           betatof_a(1) = c2/(cos(theta13)*(ds-c1))           betatof_a(1) = c2/(cos(theta13)*(ds-c1))
1087    
1088    C------- ToF Mask - S11 - S31
1089    
1090             tofmask(ch11a(tof11_i),hb11a(tof11_i)) =
1091         $        tofmask(ch11a(tof11_i),hb11a(tof11_i)) + 1
1092             tofmask(ch11b(tof11_i),hb11b(tof11_i)) =
1093         $        tofmask(ch11b(tof11_i),hb11b(tof11_i)) + 1
1094    
1095             tofmask(ch31a(tof31_i),hb31a(tof31_i)) =
1096         $        tofmask(ch31a(tof31_i),hb31a(tof31_i)) + 1
1097             tofmask(ch31b(tof31_i),hb31b(tof31_i)) =
1098         $        tofmask(ch31b(tof31_i),hb31b(tof31_i)) + 1
1099    
1100    C-------
1101          
1102        ENDIF        ENDIF
1103                
1104  C     S11 - S32  C     S11 - S32
1105        IF (tof11_i.GT.none_find.AND.tof32_i.GT.none_find) THEN  
1106           IF ((tof11_i.GT.none_find).AND.(tof32_i.GT.none_find).AND.
1107         &    (ytofpos(1).NE.101.).AND.(xtofpos(3).NE.101.)) THEN
1108           xhelp1 = tof11(1,tof11_i,itdc)+tof11(2,tof11_i,itdc)           xhelp1 = tof11(1,tof11_i,itdc)+tof11(2,tof11_i,itdc)
1109           xhelp2 = tof32(1,tof32_i,itdc)+tof32(2,tof32_i,itdc)           xhelp2 = tof32(1,tof32_i,itdc)+tof32(2,tof32_i,itdc)
1110           ds = xhelp1-xhelp2           ds = xhelp1-xhelp2
1111           ihelp=(tof11_i-1)*3+tof32_i           ihelp=(tof11_i-1)*3+tof32_i
1112           c1 = k_S11S32(1,ihelp)           c1 = k_S11S32(1,ihelp)
1113             if (iz.gt.2) c1 = c1 + k1corrA1
1114           c2 = k_S11S32(2,ihelp)           c2 = k_S11S32(2,ihelp)
1115           betatof_a(2) = c2/(cos(theta13)*(ds-c1))           betatof_a(2) = c2/(cos(theta13)*(ds-c1))
       ENDIF  
1116    
1117    C------- ToF Mask - S11 - S32
1118    
1119             tofmask(ch11a(tof11_i),hb11a(tof11_i)) =
1120         $        tofmask(ch11a(tof11_i),hb11a(tof11_i)) + 1
1121             tofmask(ch11b(tof11_i),hb11b(tof11_i)) =
1122         $        tofmask(ch11b(tof11_i),hb11b(tof11_i)) + 1
1123    
1124             tofmask(ch32a(tof32_i),hb32a(tof32_i)) =
1125         $        tofmask(ch32a(tof32_i),hb32a(tof32_i)) + 1
1126             tofmask(ch32b(tof32_i),hb32b(tof32_i)) =
1127         $        tofmask(ch32b(tof32_i),hb32b(tof32_i)) + 1
1128    
1129    C-------
1130    
1131          ENDIF
1132          
1133  C     S12 - S31  C     S12 - S31
1134        IF (tof12_i.GT.none_find.AND.tof31_i.GT.none_find) THEN  
1135           IF ((tof12_i.GT.none_find).AND.(tof31_i.GT.none_find).AND.
1136         &    (xtofpos(1).NE.101.).AND.(ytofpos(3).NE.101.)) THEN
1137           xhelp1 = tof12(1,tof12_i,itdc)+tof12(2,tof12_i,itdc)           xhelp1 = tof12(1,tof12_i,itdc)+tof12(2,tof12_i,itdc)
1138           xhelp2 = tof31(1,tof31_i,itdc)+tof31(2,tof31_i,itdc)           xhelp2 = tof31(1,tof31_i,itdc)+tof31(2,tof31_i,itdc)
1139           ds = xhelp1-xhelp2           ds = xhelp1-xhelp2
1140           ihelp=(tof12_i-1)*3+tof31_i           ihelp=(tof12_i-1)*3+tof31_i
1141           c1 = k_S12S31(1,ihelp)           c1 = k_S12S31(1,ihelp)
1142             if (iz.gt.2) c1 = c1 + k1corrA1
1143           c2 = k_S12S31(2,ihelp)           c2 = k_S12S31(2,ihelp)
1144           betatof_a(3) = c2/(cos(theta13)*(ds-c1))           betatof_a(3) = c2/(cos(theta13)*(ds-c1))
       ENDIF  
1145    
1146    C------- ToF Mask - S12 - S31
1147    
1148             tofmask(ch12a(tof12_i),hb12a(tof12_i)) =
1149         $        tofmask(ch12a(tof12_i),hb12a(tof12_i)) + 1
1150             tofmask(ch12b(tof12_i),hb12b(tof12_i)) =
1151         $        tofmask(ch12b(tof12_i),hb12b(tof12_i)) + 1
1152    
1153             tofmask(ch31a(tof31_i),hb31a(tof31_i)) =
1154         $        tofmask(ch31a(tof31_i),hb31a(tof31_i)) + 1
1155             tofmask(ch31b(tof31_i),hb31b(tof31_i)) =
1156         $        tofmask(ch31b(tof31_i),hb31b(tof31_i)) + 1
1157    
1158    C-------
1159    
1160          ENDIF
1161                
1162  C     S12 - S32  C     S12 - S32
1163        IF (tof12_i.GT.none_find.AND.tof32_i.GT.none_find) THEN  
1164           IF ((tof12_i.GT.none_find).AND.(tof32_i.GT.none_find).AND.
1165         &    (xtofpos(1).NE.101.).AND.(xtofpos(3).NE.101.)) THEN
1166           xhelp1 = tof12(1,tof12_i,itdc)+tof12(2,tof12_i,itdc)           xhelp1 = tof12(1,tof12_i,itdc)+tof12(2,tof12_i,itdc)
1167           xhelp2 = tof32(1,tof32_i,itdc)+tof32(2,tof32_i,itdc)           xhelp2 = tof32(1,tof32_i,itdc)+tof32(2,tof32_i,itdc)
1168           ds = xhelp1-xhelp2           ds = xhelp1-xhelp2
1169           ihelp=(tof12_i-1)*3+tof32_i           ihelp=(tof12_i-1)*3+tof32_i
1170           c1 = k_S12S32(1,ihelp)           c1 = k_S12S32(1,ihelp)
1171             if (iz.gt.2) c1 = c1 + k1corrA1
1172           c2 = k_S12S32(2,ihelp)           c2 = k_S12S32(2,ihelp)
1173           betatof_a(4) = c2/(cos(theta13)*(ds-c1))           betatof_a(4) = c2/(cos(theta13)*(ds-c1))
1174    
1175    C------- ToF Mask - S12 - S32
1176    
1177             tofmask(ch12a(tof12_i),hb12a(tof12_i)) =
1178         $        tofmask(ch12a(tof12_i),hb12a(tof12_i)) + 1
1179             tofmask(ch12b(tof12_i),hb12b(tof12_i)) =
1180         $        tofmask(ch12b(tof12_i),hb12b(tof12_i)) + 1
1181    
1182             tofmask(ch32a(tof32_i),hb32a(tof32_i)) =
1183         $        tofmask(ch32a(tof32_i),hb32a(tof32_i)) + 1
1184             tofmask(ch32b(tof32_i),hb32b(tof32_i)) =
1185         $        tofmask(ch32b(tof32_i),hb32b(tof32_i)) + 1
1186    
1187    C-------
1188    
1189        ENDIF        ENDIF
1190    
1191  C     S21 - S31  C     S21 - S31
1192        IF (tof21_i.GT.none_find.AND.tof31_i.GT.none_find) THEN  
1193           IF ((tof21_i.GT.none_find).AND.(tof31_i.GT.none_find).AND.
1194         &    (xtofpos(2).NE.101.).AND.(ytofpos(3).NE.101.)) THEN
1195           xhelp1 = tof21(1,tof21_i,itdc)+tof21(2,tof21_i,itdc)           xhelp1 = tof21(1,tof21_i,itdc)+tof21(2,tof21_i,itdc)
1196           xhelp2 = tof31(1,tof31_i,itdc)+tof31(2,tof31_i,itdc)           xhelp2 = tof31(1,tof31_i,itdc)+tof31(2,tof31_i,itdc)
1197           ds = xhelp1-xhelp2           ds = xhelp1-xhelp2
1198           ihelp=(tof21_i-1)*3+tof31_i           ihelp=(tof21_i-1)*3+tof31_i
1199           c1 = k_S21S31(1,ihelp)           c1 = k_S21S31(1,ihelp)
1200             if (iz.gt.2) c1 = c1 + k1corrB1
1201           c2 = k_S21S31(2,ihelp)           c2 = k_S21S31(2,ihelp)
1202           betatof_a(5) = c2/(cos(theta23)*(ds-c1))           betatof_a(5) = c2/(cos(theta13)*(ds-c1))
1203    
1204    C------- ToF Mask - S21 - S31
1205    
1206             tofmask(ch21a(tof21_i),hb21a(tof21_i)) =
1207         $        tofmask(ch21a(tof21_i),hb21a(tof21_i)) + 1
1208             tofmask(ch21b(tof21_i),hb21b(tof21_i)) =
1209         $        tofmask(ch21b(tof21_i),hb21b(tof21_i)) + 1
1210    
1211             tofmask(ch31a(tof31_i),hb31a(tof31_i)) =
1212         $        tofmask(ch31a(tof31_i),hb31a(tof31_i)) + 1
1213             tofmask(ch31b(tof31_i),hb31b(tof31_i)) =
1214         $        tofmask(ch31b(tof31_i),hb31b(tof31_i)) + 1
1215    
1216    C-------
1217    
1218        ENDIF        ENDIF
1219    
1220  C     S21 - S32  C     S21 - S32
1221        IF (tof21_i.GT.none_find.AND.tof32_i.GT.none_find) THEN  
1222           IF ((tof21_i.GT.none_find).AND.(tof32_i.GT.none_find).AND.
1223         &    (xtofpos(2).NE.101.).AND.(xtofpos(3).NE.101.)) THEN
1224           xhelp1 = tof21(1,tof21_i,itdc)+tof21(2,tof21_i,itdc)           xhelp1 = tof21(1,tof21_i,itdc)+tof21(2,tof21_i,itdc)
1225           xhelp2 = tof32(1,tof32_i,itdc)+tof32(2,tof32_i,itdc)           xhelp2 = tof32(1,tof32_i,itdc)+tof32(2,tof32_i,itdc)
1226           ds = xhelp1-xhelp2           ds = xhelp1-xhelp2
1227           ihelp=(tof21_i-1)*3+tof32_i           ihelp=(tof21_i-1)*3+tof32_i
1228           c1 = k_S21S32(1,ihelp)           c1 = k_S21S32(1,ihelp)
1229             if (iz.gt.2) c1 = c1 + k1corrB1
1230           c2 = k_S21S32(2,ihelp)           c2 = k_S21S32(2,ihelp)
1231           betatof_a(6) = c2/(cos(theta23)*(ds-c1))           betatof_a(6) = c2/(cos(theta13)*(ds-c1))
1232                      
1233    C------- ToF Mask - S21 - S32
1234    
1235             tofmask(ch21a(tof21_i),hb21a(tof21_i)) =
1236         $        tofmask(ch21a(tof21_i),hb21a(tof21_i)) + 1
1237             tofmask(ch21b(tof21_i),hb21b(tof21_i)) =
1238         $        tofmask(ch21b(tof21_i),hb21b(tof21_i)) + 1
1239    
1240             tofmask(ch32a(tof32_i),hb32a(tof32_i)) =
1241         $        tofmask(ch32a(tof32_i),hb32a(tof32_i)) + 1
1242             tofmask(ch32b(tof32_i),hb32b(tof32_i)) =
1243         $        tofmask(ch32b(tof32_i),hb32b(tof32_i)) + 1
1244    
1245    C-------
1246    
1247        ENDIF        ENDIF
1248    
1249  C     S22 - S31  C     S22 - S31
1250        IF (tof22_i.GT.none_find.AND.tof31_i.GT.none_find) THEN  
1251           IF ((tof22_i.GT.none_find).AND.(tof31_i.GT.none_find).AND.
1252         &    (ytofpos(2).NE.101.).AND.(ytofpos(3).NE.101.)) THEN
1253           xhelp1 = tof22(1,tof22_i,itdc)+tof22(2,tof22_i,itdc)           xhelp1 = tof22(1,tof22_i,itdc)+tof22(2,tof22_i,itdc)
1254           xhelp2 = tof31(1,tof31_i,itdc)+tof31(2,tof31_i,itdc)           xhelp2 = tof31(1,tof31_i,itdc)+tof31(2,tof31_i,itdc)
1255           ds = xhelp1-xhelp2           ds = xhelp1-xhelp2
1256           ihelp=(tof22_i-1)*3+tof31_i           ihelp=(tof22_i-1)*3+tof31_i
1257           c1 = k_S22S31(1,ihelp)           c1 = k_S22S31(1,ihelp)
1258             if (iz.gt.2) c1 = c1 + k1corrB1
1259           c2 = k_S22S31(2,ihelp)           c2 = k_S22S31(2,ihelp)
1260           betatof_a(7) = c2/(cos(theta13)*(ds-c1))           betatof_a(7) = c2/(cos(theta13)*(ds-c1))
1261    
1262    C------- ToF Mask - S22 - S31
1263    
1264             tofmask(ch22a(tof22_i),hb22a(tof22_i)) =
1265         $        tofmask(ch22a(tof22_i),hb22a(tof22_i)) + 1
1266             tofmask(ch22b(tof22_i),hb22b(tof22_i)) =
1267         $        tofmask(ch22b(tof22_i),hb22b(tof22_i)) + 1
1268    
1269             tofmask(ch31a(tof31_i),hb31a(tof31_i)) =
1270         $        tofmask(ch31a(tof31_i),hb31a(tof31_i)) + 1
1271             tofmask(ch31b(tof31_i),hb31b(tof31_i)) =
1272         $        tofmask(ch31b(tof31_i),hb31b(tof31_i)) + 1
1273    
1274    C-------          
1275    
1276        ENDIF        ENDIF
1277                
1278  C     S22 - S32  C     S22 - S32
1279        IF (tof22_i.GT.none_find.AND.tof32_i.GT.none_find) THEN  
1280           IF ((tof22_i.GT.none_find).AND.(tof32_i.GT.none_find).AND.
1281         &    (ytofpos(2).NE.101.).AND.(xtofpos(3).NE.101.)) THEN
1282           xhelp1 = tof22(1,tof22_i,itdc)+tof22(2,tof22_i,itdc)           xhelp1 = tof22(1,tof22_i,itdc)+tof22(2,tof22_i,itdc)
1283           xhelp2 = tof32(1,tof32_i,itdc)+tof32(2,tof32_i,itdc)           xhelp2 = tof32(1,tof32_i,itdc)+tof32(2,tof32_i,itdc)
1284           ds = xhelp1-xhelp2           ds = xhelp1-xhelp2
1285           ihelp=(tof22_i-1)*3+tof32_i           ihelp=(tof22_i-1)*3+tof32_i
1286           c1 = k_S22S32(1,ihelp)           c1 = k_S22S32(1,ihelp)
1287             if (iz.gt.2) c1 = c1 + k1corrB1
1288           c2 = k_S22S32(2,ihelp)           c2 = k_S22S32(2,ihelp)
1289           betatof_a(8) = c2/(cos(theta13)*(ds-c1))           betatof_a(8) = c2/(cos(theta13)*(ds-c1))
1290    
1291    C------- ToF Mask - S22 - S32
1292    
1293             tofmask(ch22a(tof22_i),hb22a(tof22_i)) =
1294         $        tofmask(ch22a(tof22_i),hb22a(tof22_i)) + 1
1295             tofmask(ch22b(tof22_i),hb22b(tof22_i)) =
1296         $        tofmask(ch22b(tof22_i),hb22b(tof22_i)) + 1
1297    
1298             tofmask(ch32a(tof32_i),hb32a(tof32_i)) =
1299         $        tofmask(ch32a(tof32_i),hb32a(tof32_i)) + 1
1300             tofmask(ch32b(tof32_i),hb32b(tof32_i)) =
1301         $        tofmask(ch32b(tof32_i),hb32b(tof32_i)) + 1
1302    
1303    C-------  
1304    
1305        ENDIF        ENDIF
1306    
1307  C     S11 - S21  C     S11 - S21
1308        IF (tof11_i.GT.none_find.AND.tof21_i.GT.none_find) THEN  
1309           IF ((tof11_i.GT.none_find).AND.(tof21_i.GT.none_find).AND.
1310         &    (ytofpos(1).NE.101.).AND.(xtofpos(2).NE.101.)) THEN
1311           xhelp1 = tof11(1,tof11_i,itdc)+tof11(2,tof11_i,itdc)           xhelp1 = tof11(1,tof11_i,itdc)+tof11(2,tof11_i,itdc)
1312           xhelp2 = tof21(1,tof21_i,itdc)+tof21(2,tof21_i,itdc)           xhelp2 = tof21(1,tof21_i,itdc)+tof21(2,tof21_i,itdc)
1313           ds = xhelp1-xhelp2           ds = xhelp1-xhelp2
1314           ihelp=(tof11_i-1)*2+tof21_i           ihelp=(tof11_i-1)*2+tof21_i
1315           c1 = k_S11S21(1,ihelp)           c1 = k_S11S21(1,ihelp)
1316             if (iz.gt.2) c1 = c1 + k1corrC1
1317           c2 = k_S11S21(2,ihelp)           c2 = k_S11S21(2,ihelp)
1318           betatof_a(9) = c2/(cos(theta13)*(ds-c1))           betatof_a(9) = c2/(cos(theta13)*(ds-c1))
1319    
1320    C------- ToF Mask - S11 - S21
1321    
1322             tofmask(ch11a(tof11_i),hb11a(tof11_i)) =
1323         $        tofmask(ch11a(tof11_i),hb11a(tof11_i)) + 1
1324             tofmask(ch11b(tof11_i),hb11b(tof11_i)) =
1325         $        tofmask(ch11b(tof11_i),hb11b(tof11_i)) + 1
1326    
1327             tofmask(ch21a(tof21_i),hb21a(tof21_i)) =
1328         $        tofmask(ch21a(tof21_i),hb21a(tof21_i)) + 1
1329             tofmask(ch21b(tof21_i),hb21b(tof21_i)) =
1330         $        tofmask(ch21b(tof21_i),hb21b(tof21_i)) + 1
1331    
1332    C-------  
1333    
1334        ENDIF        ENDIF
1335                
1336  C     S11 - S22  C     S11 - S22
1337        IF (tof11_i.GT.none_find.AND.tof22_i.GT.none_find) THEN  
1338           IF ((tof11_i.GT.none_find).AND.(tof22_i.GT.none_find).AND.
1339         &    (ytofpos(1).NE.101.).AND.(ytofpos(2).NE.101.)) THEN
1340           xhelp1 = tof11(1,tof11_i,itdc)+tof11(2,tof11_i,itdc)           xhelp1 = tof11(1,tof11_i,itdc)+tof11(2,tof11_i,itdc)
1341           xhelp2 = tof22(1,tof22_i,itdc)+tof22(2,tof22_i,itdc)           xhelp2 = tof22(1,tof22_i,itdc)+tof22(2,tof22_i,itdc)
1342           ds = xhelp1-xhelp2           ds = xhelp1-xhelp2
1343           ihelp=(tof11_i-1)*2+tof22_i           ihelp=(tof11_i-1)*2+tof22_i
1344           c1 = k_S11S22(1,ihelp)           c1 = k_S11S22(1,ihelp)
1345             if (iz.gt.2) c1 = c1 + k1corrC1
1346           c2 = k_S11S22(2,ihelp)           c2 = k_S11S22(2,ihelp)
1347           betatof_a(10) = c2/(cos(theta13)*(ds-c1))           betatof_a(10) = c2/(cos(theta13)*(ds-c1))
1348    
1349    C------- ToF Mask - S11 - S22
1350    
1351             tofmask(ch11a(tof11_i),hb11a(tof11_i)) =
1352         $        tofmask(ch11a(tof11_i),hb11a(tof11_i)) + 1
1353             tofmask(ch11b(tof11_i),hb11b(tof11_i)) =
1354         $        tofmask(ch11b(tof11_i),hb11b(tof11_i)) + 1
1355    
1356             tofmask(ch22a(tof22_i),hb22a(tof22_i)) =
1357         $        tofmask(ch22a(tof22_i),hb22a(tof22_i)) + 1
1358             tofmask(ch22b(tof22_i),hb22b(tof22_i)) =
1359         $        tofmask(ch22b(tof22_i),hb22b(tof22_i)) + 1
1360    
1361    C-------  
1362    
1363        ENDIF        ENDIF
1364    
1365  C     S12 - S21  C     S12 - S21
1366        IF (tof12_i.GT.none_find.AND.tof21_i.GT.none_find) THEN  
1367           IF ((tof12_i.GT.none_find).AND.(tof21_i.GT.none_find).AND.
1368         &    (xtofpos(1).NE.101.).AND.(xtofpos(2).NE.101.)) THEN
1369           xhelp1 = tof12(1,tof12_i,itdc)+tof12(2,tof12_i,itdc)           xhelp1 = tof12(1,tof12_i,itdc)+tof12(2,tof12_i,itdc)
1370           xhelp2 = tof21(1,tof21_i,itdc)+tof21(2,tof21_i,itdc)           xhelp2 = tof21(1,tof21_i,itdc)+tof21(2,tof21_i,itdc)
1371           ds = xhelp1-xhelp2           ds = xhelp1-xhelp2
1372           ihelp=(tof12_i-1)*2+tof21_i           ihelp=(tof12_i-1)*2+tof21_i
1373           c1 = k_S12S21(1,ihelp)           c1 = k_S12S21(1,ihelp)
1374             if (iz.gt.2) c1 = c1 + k1corrC1
1375           c2 = k_S12S21(2,ihelp)           c2 = k_S12S21(2,ihelp)
1376           betatof_a(11) = c2/(cos(theta13)*(ds-c1))           betatof_a(11) = c2/(cos(theta13)*(ds-c1))
1377    
1378    C------- ToF Mask - S12 - S21
1379    
1380             tofmask(ch12a(tof12_i),hb12a(tof12_i)) =
1381         $        tofmask(ch12a(tof12_i),hb12a(tof12_i)) + 1
1382             tofmask(ch12b(tof12_i),hb12b(tof12_i)) =
1383         $        tofmask(ch12b(tof12_i),hb12b(tof12_i)) + 1
1384    
1385             tofmask(ch21a(tof21_i),hb21a(tof21_i)) =
1386         $        tofmask(ch21a(tof21_i),hb21a(tof21_i)) + 1
1387             tofmask(ch21b(tof21_i),hb21b(tof21_i)) =
1388         $        tofmask(ch21b(tof21_i),hb21b(tof21_i)) + 1
1389    
1390    C-------  
1391    
1392        ENDIF        ENDIF
1393    
1394  C     S12 - S22  C     S12 - S22
1395        IF (tof12_i.GT.none_find.AND.tof22_i.GT.none_find) THEN  
1396           IF ((tof12_i.GT.none_find).AND.(tof22_i.GT.none_find).AND.
1397         &    (xtofpos(1).NE.101.).AND.(ytofpos(2).NE.101.)) THEN
1398           xhelp1 = tof12(1,tof12_i,itdc)+tof12(2,tof12_i,itdc)           xhelp1 = tof12(1,tof12_i,itdc)+tof12(2,tof12_i,itdc)
1399           xhelp2 = tof22(1,tof22_i,itdc)+tof22(2,tof22_i,itdc)           xhelp2 = tof22(1,tof22_i,itdc)+tof22(2,tof22_i,itdc)
1400           ds = xhelp1-xhelp2           ds = xhelp1-xhelp2
1401           ihelp=(tof12_i-1)*2+tof22_i           ihelp=(tof12_i-1)*2+tof22_i
1402           c1 = k_S12S22(1,ihelp)           c1 = k_S12S22(1,ihelp)
1403             if (iz.gt.2) c1 = c1 + k1corrC1
1404           c2 = k_S12S22(2,ihelp)           c2 = k_S12S22(2,ihelp)
1405           betatof_a(12) = c2/(cos(theta13)*(ds-c1))           betatof_a(12) = c2/(cos(theta13)*(ds-c1))
1406    
1407    C------- ToF Mask - S12 - S22
1408    
1409             tofmask(ch12a(tof12_i),hb12a(tof12_i)) =
1410         $        tofmask(ch12a(tof12_i),hb12a(tof12_i)) + 1
1411             tofmask(ch12b(tof12_i),hb12b(tof12_i)) =
1412         $        tofmask(ch12b(tof12_i),hb12b(tof12_i)) + 1
1413    
1414             tofmask(ch22a(tof22_i),hb22a(tof22_i)) =
1415         $        tofmask(ch22a(tof22_i),hb22a(tof22_i)) + 1
1416             tofmask(ch22b(tof22_i),hb22b(tof22_i)) =
1417         $        tofmask(ch22b(tof22_i),hb22b(tof22_i)) + 1
1418    
1419    C-------  
1420    
1421        ENDIF        ENDIF
         
 C-------  
1422    
1423        icount=0  C---------------------------------------------------------
1424        sw=0.  C
1425        sxw=0.  C      icount=0
1426        beta_mean=100.  C      sw=0.
1427    C      sxw=0.
1428    C      beta_mean=100.
1429    C
1430    C      do i=1,12
1431    C         if ((betatof_a(i).gt.-1.5).and.(betatof_a(i).lt.1.5)) then
1432    C            icount= icount+1
1433    C            if (i.le.4) w_i=1./(0.13**2.)
1434    C            if ((i.ge.5).and.(i.le.8)) w_i=1./(0.16**2.)
1435    C         if (i.ge.9) w_i=1./(0.25**2.)     ! to be checked
1436    C            sxw=sxw + betatof_a(i)*w_i
1437    C            sw =sw + w_i
1438    C         endif
1439    C      enddo
1440    C      
1441    C      if (icount.gt.0) beta_mean=sxw/sw
1442    C      betatof_a(13) = beta_mean
1443    C
1444    
1445    C--------  New mean beta  calculation  -----------------------
1446    
1447        do i=1,12          do i=1,12
1448           if ((betatof_a(i).gt.-1.5).and.(betatof_a(i).lt.1.5)) then           btemp(i) = betatof_a(i)
1449              icount= icount+1          enddo
1450              if (i.le.4) w_i=1./(0.13**2.)  
1451              if ((i.ge.5).and.(i.le.8)) w_i=1./(0.16**2.)          betatof_a(13)=newbeta(1,btemp,hitvec,10.,10.,20.)
1452              if (i.ge.9) w_i=1./(0.25**2.)     ! to be checked  
1453              sxw=sxw + betatof_a(i)*w_i  C--------------------------------------------------------------
1454              sw =sw + w_i  C      write(*,*) betatof_a
1455           endif  c      write(*,*) xtofpos
1456        enddo  c      write(*,*) ytofpos
1457    c      write(*,*)'tofl2com beta', betatof_a
1458    C      write(*,*) adcflagtof
1459    c      write(*,*) 'tofl2com'
1460    c      write(*,*) xtofpos
1461    c      write(*,*) ytofpos
1462    c      write(*,*) xtr_tof
1463    c      write(*,*) ytr_tof
1464                
1465        if (icount.gt.0) beta_mean=sxw/sw  c 100  continue
1466        betatof_a(13) = beta_mean          continue
 c      write(*,*) '------------- end  tofl2com ----------'  
   
  100  continue  
1467    
1468  C  C
1469        RETURN        RETURN
1470        END        END
1471    
1472    
1473    C------------------------------------------------------------------
1474    C------------------------------------------------------------------
1475    
1476           function atten(is,ilay,ipad,x)
1477           include  'input_tof.txt'
1478           real atten
1479           real x
1480           real xmin,xmax
1481           integer  ilay,ipad
1482    
1483    *  S11 8 paddles  33.0 x 5.1 cm
1484    *  S12 6 paddles  40.8 x 5.5 cm
1485    *  S21 2 paddles  18.0 x 7.5 cm
1486    *  S22 2 paddles  15.0 x 9.0 cm
1487    *  S31 3 paddles  15.0 x 6.0 cm
1488    *  S32 3 paddles  18.0 x 5.0 cm
1489    
1490    
1491    c       if (ilay.eq.11) write(*,*) 'start ',ipad,is,adcx11(is,ipad,1),
1492    c     &  adcx11(is,ipad,2),adcx11(is,ipad,3),adcx11(is,ipad,4)
1493    c       if (ilay.eq.12) write(*,*) 'start ',ipad,is,adcx12(is,ipad,1),
1494    c     &  adcx12(is,ipad,2),adcx12(is,ipad,3),adcx12(is,ipad,4)
1495    
1496    
1497           if (ilay.eq.11)  xmin=-33.0/2.
1498           if (ilay.eq.11)  xmax= 33.0/2.
1499           if (ilay.eq.12)  xmin=-40.8/2.
1500           if (ilay.eq.12)  xmax= 40.8/2.
1501    
1502           if (ilay.eq.21)  xmin=-18.0/2.
1503           if (ilay.eq.21)  xmax= 18.0/2.
1504           if (ilay.eq.22)  xmin=-15.0/2.
1505           if (ilay.eq.22)  xmax= 15.0/2.
1506    
1507           if (ilay.eq.31)  xmin=-15.0/2.
1508           if (ilay.eq.31)  xmax= 15.0/2.
1509           if (ilay.eq.32)  xmin=-18.0/2.
1510           if (ilay.eq.32)  xmax= 18.0/2.
1511    
1512           if (x .lt. xmin) x=xmin
1513           if (x .gt. xmax) x=xmax
1514    
1515    
1516           if (ilay.eq.11) atten=
1517         &    adcx11(is,ipad,1)*exp(x*adcx11(is,ipad,2))
1518         &  + adcx11(is,ipad,3)*exp(x*adcx11(is,ipad,4))
1519    
1520           if (ilay.eq.12) atten=
1521         &    adcx12(is,ipad,1)*exp(x*adcx12(is,ipad,2))
1522         &  + adcx12(is,ipad,3)*exp(x*adcx12(is,ipad,4))
1523    
1524           if (ilay.eq.21) atten=
1525         &    adcx21(is,ipad,1)*exp(x*adcx21(is,ipad,2))
1526         &  + adcx21(is,ipad,3)*exp(x*adcx21(is,ipad,4))
1527    
1528           if (ilay.eq.22) atten=
1529         &    adcx22(is,ipad,1)*exp(x*adcx22(is,ipad,2))
1530         &  + adcx22(is,ipad,3)*exp(x*adcx22(is,ipad,4))
1531    
1532           if (ilay.eq.31) atten=
1533         &    adcx31(is,ipad,1)*exp(x*adcx31(is,ipad,2))
1534         &  + adcx31(is,ipad,3)*exp(x*adcx31(is,ipad,4))
1535    
1536           if (ilay.eq.32) atten=
1537         &    adcx32(is,ipad,1)*exp(x*adcx32(is,ipad,2))
1538         &  + adcx32(is,ipad,3)*exp(x*adcx32(is,ipad,4))
1539    
1540            if (atten.gt.10000) atten=10000.
1541    
1542           end
1543    
1544    C------------------------------------------------------------------
1545    C------------------------------------------------------------------
1546    
1547           function pc_adc(ix)
1548           include  'input_tof.txt'
1549           real pc_adc
1550           integer ix
1551    
1552           pc_adc=28.0407 + 0.628929*ix
1553         &   - 5.80901e-05*ix*ix + 3.14092e-08*ix*ix*ix
1554    c       write(*,*) ix,pc_adc
1555           end
1556    
1557    C------------------------------------------------------------------
1558    C------------------------------------------------------------------
1559    
1560            function check_charge(theta,hitvec)
1561    
1562            include  'input_tof.txt'
1563            include  'tofcomm.txt'
1564    
1565            real check_charge  
1566            integer hitvec(6)  
1567            REAL CHARGE, theta
1568    
1569    C  upper and lower limits  for the helium selection
1570            REAL A_l(24),A_h(24)
1571            DATA A_l /200,190,300,210,220,200,210,60,60,120,220,
1572         &  120,160,50,300,200,120,250,350,300,350,250,280,300/
1573            DATA A_h /550,490,800,600,650,600,600,260,200,380,
1574         &  620,380,550,200,850,560,400,750,900,800,880,800,750,800/
1575    
1576    C The k1 constants for the beta calculation, only for S1-S3
1577    C k2 constant is taken to be the standard 2D/c
1578            REAL k1(84)
1579            DATA k1 /50,59.3296,28.4328,-26.0818,5.91253,-19.588,
1580         &   -9.26316,24.7544,2.32465,-50.5058,-15.3195,-39.1443,
1581         &   -91.2546,-58.6243,-84.5641,-63.1516,-32.2091,-58.3358,
1582         &   13.8084,45.5322,33.2416,-11.5313,51.3271,75,-14.1141,
1583         &   42.8466,15.1794,-63.6672,-6.07739,-32.164,-41.771,10.5274,
1584         &   -9.46096,-81.7404,-28.783,-52.7167,-127.394,-69.6166,
1585         &   -93.4655,-98.9543,-42.863,-67.8244,-19.3238,31.1221,8.7319,
1586         &   -43.1627,5.55573,-14.4078,-83.4466,-47.4647,-77.8379,
1587         &   -108.222,-75.986,-101.297,-96.0205,-63.1881,-90.1372,
1588         &   -22.7347,8.31409,-19.6912,-7.49008,23.6979,-1.66677,
1589         &   1.81556,34.4668,6.23693,-100,-59.5861,-90.9159,-141.639,
1590         &   -89.2521,-112.881,-130.199,-77.0357,-98.4632,-60.2086,
1591         &   -4.82097,-29.3705,-43.6469,10.5884,-9.31304,-35.3329,
1592         &   25.2514,25.6/
1593    
1594    
1595    
1596            REAL zin(6)
1597            DATA zin /53.74, 53.04, 23.94, 23.44, -23.49, -24.34/
1598    
1599            REAL  c1,c2,xhelp,xhelp1,xhelp2,ds,dist,F
1600            REAL  sw,sxw,beta_mean_tof,w_i
1601            INTEGER  ihelp
1602            INTEGER ipmt(4)
1603            REAL time(4),beta1(4)
1604    
1605            REAL  adca(48), tdca(48)
1606    
1607            REAL  a1,a2
1608            INTEGER jj
1609    
1610    c get rid of warnings EMILIANO
1611            i = 0
1612            slope = 0
1613            offset = 0
1614            none_find = 0
1615            none_ev = 0
1616            adc_ev = 0
1617            tdc_ev = 0
1618            iadc = 0
1619            itdc = 0
1620            right = 0
1621            left = 0
1622            tof12_y(1) = tof12_y(1)
1623            tof11_x(1) = tof11_x(1)
1624            tof21_y(1) = tof21_y(1)
1625            tof22_x(1) = tof22_x(1)
1626            tof32_y(1) = tof32_y(1)
1627            tof31_x(1) = tof31_x(1)
1628    c get rid of warnings
1629    
1630    C-----------------------------------------------------------
1631    C--- get data
1632    C-----------------------------------------------------------
1633    
1634             do j=1,8
1635             ih = 1 + 0 +((j-1)*2)
1636             adca(ih)   = adc(ch11a(j),hb11a(j))
1637             adca(ih+1) = adc(ch11b(j),hb11b(j))
1638             tdca(ih)   = tdc(ch11a(j),hb11a(j))
1639             tdca(ih+1) = tdc(ch11b(j),hb11b(j))
1640             enddo
1641    
1642             do j=1,6
1643             ih = 1 + 16+((j-1)*2)
1644             adca(ih)   = adc(ch12a(j),hb12a(j))
1645             adca(ih+1) = adc(ch12b(j),hb12b(j))
1646             tdca(ih)   = tdc(ch12a(j),hb12a(j))
1647             tdca(ih+1) = tdc(ch12b(j),hb12b(j))
1648             enddo
1649    
1650             do j=1,2
1651             ih = 1 + 28+((j-1)*2)
1652             adca(ih)   = adc(ch21a(j),hb21a(j))
1653             adca(ih+1) = adc(ch21b(j),hb21b(j))
1654             tdca(ih)   = tdc(ch21a(j),hb21a(j))
1655             tdca(ih+1) = tdc(ch21b(j),hb21b(j))
1656             enddo
1657    
1658             do j=1,2
1659             ih = 1 + 32+((j-1)*2)
1660             adca(ih)   = adc(ch22a(j),hb22a(j))
1661             adca(ih+1) = adc(ch22b(j),hb22b(j))
1662             tdca(ih)   = tdc(ch22a(j),hb22a(j))
1663             tdca(ih+1) = tdc(ch22b(j),hb22b(j))
1664             enddo
1665    
1666             do j=1,3
1667             ih = 1 + 36+((j-1)*2)
1668             adca(ih)   = adc(ch31a(j),hb31a(j))
1669             adca(ih+1) = adc(ch31b(j),hb31b(j))
1670             tdca(ih)   = tdc(ch31a(j),hb31a(j))
1671             tdca(ih+1) = tdc(ch31b(j),hb31b(j))
1672             enddo
1673    
1674             do j=1,3
1675             ih = 1 + 42+((j-1)*2)
1676             adca(ih)   = adc(ch32a(j),hb32a(j))
1677             adca(ih+1) = adc(ch32b(j),hb32b(j))
1678             tdca(ih)   = tdc(ch32a(j),hb32a(j))
1679             tdca(ih+1) = tdc(ch32b(j),hb32b(j))
1680             enddo
1681    
1682    
1683    c         write(*,*) adca
1684    c         write(*,*) tdca
1685    
1686    
1687    C============   calculate beta and select charge > Z=1   ===============
1688    
1689            ICHARGE=1
1690    
1691    C find hitted paddle by looking for ADC values on both sides
1692    C since we looking for Z>1 this gives decent results
1693    
1694            tof11_i = hitvec(1)-1
1695            tof12_i = hitvec(2)-1
1696            tof21_i = hitvec(3)-1
1697            tof22_i = hitvec(4)-1
1698            tof31_i = hitvec(5)-1
1699            tof32_i = hitvec(6)-1
1700    
1701    c        write(*,*) ' in charge check'
1702    c        write(*,*) theta,tof11_i,tof12_i,tof21_i,tof22_i,tof31_i,tof32_i
1703    
1704    C----------------------------------------------------------------
1705    
1706            beta_help=100.
1707            beta_mean_tof=100.
1708    
1709            do jj=1,4
1710              beta1(jj) = 100.
1711            enddo
1712    
1713    C----------------------------------------------------------------
1714    C---------  S1 - S3 ---------------------------------------------
1715    C----------------------------------------------------------------
1716    
1717    C---------  S11 - S31 -------------------------------------------
1718    
1719            if ((tof11_i.gt.-1).and.(tof31_i.gt.-1)) then
1720    
1721            dist = zin(1) - zin(5)
1722            c2 = (2.*0.01*dist)/(3.E08*50.E-12)
1723            F = 1./cos(theta)
1724    
1725            ipmt(1)   = (tof11_i)*2+1
1726            ipmt(2)   = (tof11_i)*2+2
1727            ipmt(3)   = 36+(tof31_i)*2+1
1728            ipmt(4)   = 36+(tof31_i)*2+2
1729    
1730    c        write(*,*) ipmt
1731    
1732            do jj=1,4
1733               time(jj) = tdca(ipmt(jj))
1734            enddo
1735    
1736    c        write(*,*) time
1737    
1738            if ((time(1).lt.4095).and.(time(2).lt.4095).and.
1739         &     (time(3).lt.4095).and.(time(4).lt.4095)) then
1740             xhelp1 = time(1) + time(2)
1741             xhelp2 = time(3) + time(4)
1742             ds = xhelp1-xhelp2
1743             ihelp=0+(tof11_i)*3+tof31_i
1744             c1 = k1(ihelp+1)
1745             beta1(1) = c2*F/(ds-c1);
1746                                                     endif
1747    c         write(*,*) beta1(1)
1748             endif  ! tof_....
1749    
1750    
1751    C---------  S11 - S32 -------------------------------------------
1752    
1753            if ((tof11_i.gt.-1).and.(tof32_i.gt.-1)) then
1754    
1755            dist = zin(1) - zin(6)
1756            c2 = (2.*0.01*dist)/(3.E08*50.E-12)
1757            F = 1./cos(theta)
1758    
1759            ipmt(1)   = (tof11_i)*2+1
1760            ipmt(2)   = (tof11_i)*2+2
1761            ipmt(3)   = 42+(tof32_i)*2+1
1762            ipmt(4)   = 42+(tof32_i)*2+2
1763    
1764            do jj=1,4
1765               time(jj) = tdca(ipmt(jj))
1766            enddo
1767    
1768            if ((time(1).lt.4095).and.(time(2).lt.4095).and.
1769         &     (time(3).lt.4095).and.(time(4).lt.4095)) then
1770             xhelp1 = time(1) + time(2)
1771             xhelp2 = time(3) + time(4)
1772             ds = xhelp1-xhelp2
1773             ihelp=24+(tof11_i)*3+tof32_i
1774             c1 = k1(ihelp+1)
1775             beta1(2) = c2*F/(ds-c1);
1776                                                     endif
1777             endif  ! tof_....
1778    
1779    
1780    C---------  S12 - S31 -------------------------------------------
1781    
1782            if ((tof12_i.gt.-1).and.(tof31_i.gt.-1)) then
1783    
1784            dist = zin(2) - zin(5)
1785            c2 = (2.*0.01*dist)/(3.E08*50.E-12)
1786            F = 1./cos(theta)
1787    
1788            ipmt(1)   = 16+(tof12_i)*2+1
1789            ipmt(2)   = 16+(tof12_i)*2+2
1790            ipmt(3)   = 36+(tof31_i)*2+1
1791            ipmt(4)   = 36+(tof31_i)*2+2
1792    
1793            do jj=1,4
1794               time(jj) = tdca(ipmt(jj))
1795            enddo
1796    
1797            if ((time(1).lt.4095).and.(time(2).lt.4095).and.
1798         &     (time(3).lt.4095).and.(time(4).lt.4095)) then
1799             xhelp1 = time(1) + time(2)
1800             xhelp2 = time(3) + time(4)
1801             ds = xhelp1-xhelp2
1802             ihelp=48+(tof12_i)*3+tof31_i
1803             c1 = k1(ihelp+1)
1804             beta1(3) = c2*F/(ds-c1);
1805                                                     endif
1806             endif  ! tof_....
1807    
1808    
1809    C---------  S12 - S32 -------------------------------------------
1810    
1811            if ((tof12_i.gt.-1).and.(tof32_i.gt.-1)) then
1812    
1813            dist = zin(2) - zin(6)
1814            c2 = (2.*0.01*dist)/(3.E08*50.E-12)
1815            F = 1./cos(theta)
1816    
1817            ipmt(1)   = 16+(tof12_i)*2+1
1818            ipmt(2)   = 16+(tof12_i)*2+2
1819            ipmt(3)   = 42+(tof32_i)*2+1
1820            ipmt(4)   = 42+(tof32_i)*2+2
1821    
1822            do jj=1,4
1823               time(jj) = tdca(ipmt(jj))
1824            enddo
1825    
1826            if ((time(1).lt.4095).and.(time(2).lt.4095).and.
1827         &     (time(3).lt.4095).and.(time(4).lt.4095)) then
1828             xhelp1 = time(1) + time(2)
1829             xhelp2 = time(3) + time(4)
1830             ds = xhelp1-xhelp2
1831             ihelp=56+(tof12_i)*3+tof32_i
1832             c1 = k1(ihelp+1)
1833             beta1(4) = c2*F/(ds-c1);
1834                                                     endif
1835    
1836             endif  ! tof_....
1837    
1838    c         write(*,*) beta1
1839      
1840    C----  calculate  beta mean, only downward going particles are interesting ----
1841    
1842             sw=0.
1843             sxw=0.
1844             beta_mean_tof=100.
1845    
1846            do jj=1,4
1847            if ((beta1(jj).gt.0.1).and.(beta1(jj).lt.2.0)) then
1848                w_i=1./(0.13*0.13)
1849                sxw=sxw + beta1(jj)*w_i
1850                sw =sw + w_i ;
1851                                                          endif
1852            enddo
1853    
1854            if (sw.gt.0) beta_mean_tof=sxw/sw;
1855    
1856    c        write(*,*) 'beta_mean_tof ',beta_mean_tof
1857    
1858            beta_help = beta_mean_tof  !  pow(beta_mean_tof,1.0) gave best results
1859    
1860    CCCCC        endif  !  if tof11_i > -1 && ...... beta calculation
1861    
1862    C-----------------------  Select charge   --------------------------
1863    
1864           charge=0
1865    
1866           if ((beta_mean_tof.gt.0.2).and.(beta_mean_tof.lt.2.0)) then
1867    
1868            icount1=0
1869            icount2=0
1870            icount3=0
1871    
1872            do jj=0,23
1873            a1 = adca(2*jj+1)
1874            a2 = adca(2*jj+2)
1875            if ((a1.lt.4095).and.(a2.lt.4095)) then
1876            a1 = adca(2*jj+1)*cos(theta)
1877            a2 = adca(2*jj+2)*cos(theta)
1878            xhelp  = 100000.
1879            xhelp1 = 100000.
1880            xhelp = sqrt(a1*a2)  ! geometric mean
1881            xhelp1 = beta_help*xhelp
1882    C if geometric mean multiplied by beta_help is inside/outside helium
1883    C limits, increase counter
1884           if (xhelp1.lt.A_l(jj+1))  icount1=icount1+1
1885           if ((xhelp1.gt.A_l(jj+1)).and.(xhelp1.lt.A_h(jj+1)))
1886         &                           icount2=icount2+1
1887           if (xhelp1.gt.A_h(jj+1))  icount3=icount3+1
1888                                                endif
1889            enddo
1890    
1891    
1892    C  if more than three paddles see the same...
1893    
1894            if (icount1 .gt. 3) charge=1
1895            if (icount2 .gt. 3) charge=2
1896            if (icount3 .gt. 3) charge=3
1897    
1898                                                        endif  ! 0.2<beta<2.0
1899    
1900    C  no beta found? Sum up geometric means of paddles and derive the mean...
1901    
1902           if (beta_mean_tof.eq.100.) then
1903    
1904           xhelp  = 0.
1905           icount = 0
1906    
1907           if (tof11_i.gt.-1) then
1908           jj=tof11_i
1909           a1 = adca(0+2*jj+1)
1910           a2 = adca(0+2*jj+2)
1911           if ((a1.lt.4095).and.(a2.lt.4095)) then
1912           a1 = a1*cos(theta)
1913           a2 = a2*cos(theta)
1914           xhelp = xhelp + sqrt(a1*a2)
1915           icount=icount+1
1916                        endif
1917                        endif
1918    
1919           if (tof12_i.gt.-1) then
1920           jj=tof12_i
1921           a1 = adca(16+2*jj+1)
1922           a2 = adca(16+2*jj+2)
1923           if ((a1.lt.4095).and.(a2.lt.4095)) then
1924           a1 = a1*cos(theta)
1925           a2 = a2*cos(theta)
1926           xhelp = xhelp + sqrt(a1*a2)
1927           icount=icount+1
1928                        endif
1929                        endif
1930    
1931           if (tof21_i.gt.-1) then
1932           jj=tof21_i
1933           a1 = adca(28+2*jj+1)
1934           a2 = adca(28+2*jj+2)
1935           if ((a1.lt.4095).and.(a2.lt.4095)) then
1936           a1 = a1*cos(theta)
1937           a2 = a2*cos(theta)
1938           xhelp = xhelp + sqrt(a1*a2)
1939           icount=icount+1
1940                        endif
1941                        endif
1942    
1943           if (tof22_i.gt.-1) then
1944           jj=tof22_i
1945           a1 = adca(32+2*jj+1)
1946           a2 = adca(32+2*jj+2)
1947           if ((a1.lt.4095).and.(a2.lt.4095)) then
1948           a1 = a1*cos(theta)
1949           a2 = a2*cos(theta)
1950           xhelp = xhelp + sqrt(a1*a2)
1951           icount=icount+1
1952                        endif
1953                        endif
1954    
1955           if (tof31_i.gt.-1) then
1956           jj=tof31_i
1957           a1 = adca(36+2*jj+1)
1958           a2 = adca(36+2*jj+2)
1959           if ((a1.lt.4095).and.(a2.lt.4095)) then
1960           a1 = a1*cos(theta)
1961           a2 = a2*cos(theta)
1962           xhelp = xhelp + sqrt(a1*a2)
1963           icount=icount+1
1964                        endif
1965                        endif
1966    
1967           if (tof32_i.gt.-1) then
1968           jj=tof32_i
1969           a1 = adca(42+2*jj+1)
1970           a2 = adca(42+2*jj+2)
1971           if ((a1.lt.4095).and.(a2.lt.4095)) then
1972           a1 = a1*cos(theta)
1973           a2 = a2*cos(theta)
1974           xhelp = xhelp + sqrt(a1*a2)
1975           icount=icount+1
1976                        endif
1977                        endif
1978    
1979    
1980           if (icount.gt.0) xhelp=xhelp/icount
1981           if ((icount.gt.2).and.(xhelp.gt.1500.)) charge=3
1982    
1983                                      endif  ! beta_mean_tof.eq.100.
1984    
1985    c        write(*,*) 'in function charge: ',charge
1986            check_charge = charge
1987    
1988    
1989            END
1990    
1991    C****************************************************************************
1992    C****************************************************************************
1993    C****************************************************************************
1994    
1995            function newbeta(iflag,b,hitvec,resmax,qualitycut,chi2cut)
1996    
1997            include  'input_tof.txt'
1998            include  'output_tof.txt'
1999            include  'tofcomm.txt'
2000    
2001            REAL newbeta
2002            REAL resmax,qualitycut,chi2cut
2003            REAL w_i(12),w_il(6),quality,res,betachi,beta_mean_inv
2004            REAL sw,sxw,b(12),beta_mean,chi2,xhelp
2005            REAL tdcfl(4,12)
2006    
2007            INTEGER iflag,icount,hitvec(6)
2008    
2009            INTEGER itop(12),ibot(12)
2010            DATA itop /1,1,2,2,3,3,4,4,1,1,2,2/
2011            DATA ibot /5,6,5,6,5,6,5,6,3,4,3,4/
2012    
2013    
2014    c get rid of warnings EMILIANO
2015            slope = 0
2016            offset = 0
2017            none_find = 0
2018            none_ev = 0
2019            adc_ev = 0
2020            tdc_ev = 0
2021            iadc = 0
2022            itdc = 0
2023            right = 0
2024            left = 0
2025            tof12_y(1) = tof12_y(1)
2026            tof11_x(1) = tof11_x(1)
2027            tof21_y(1) = tof21_y(1)
2028            tof22_x(1) = tof22_x(1)
2029            tof32_y(1) = tof32_y(1)
2030            tof31_x(1) = tof31_x(1)
2031    c get rid of warnings
2032    
2033    C====================================================================
2034    
2035            tof11_i = hitvec(1)
2036            tof12_i = hitvec(2)
2037            tof21_i = hitvec(3)
2038            tof22_i = hitvec(4)
2039            tof31_i = hitvec(5)
2040            tof32_i = hitvec(6)
2041    
2042             if (iflag.eq.1) then   ! call from tofl2com
2043             do i=1,4
2044             do j=1,12
2045              tdcfl(i,j) =  tdcflagtof(i,j)
2046             enddo
2047             enddo
2048                            endif
2049    
2050             if (iflag.eq.2) then   ! call from toftrk
2051             do i=1,4
2052             do j=1,12
2053              tdcfl(i,j) =  tdcflag(i,j)
2054             enddo
2055             enddo
2056                            endif
2057    
2058    
2059    C---  Find out ToF layers with artificial TDC values    -------------
2060    
2061            do jj=1,6
2062            w_il(jj) = 1000.
2063            enddo
2064    
2065    
2066            if (tof11_i.gt.0) then
2067            if ((tofmask(ch11a(tof11_i),hb11a(tof11_i)).gt.0).or.
2068         &   (tofmask(ch11b(tof11_i),hb11b(tof11_i)).gt.0)) then
2069            w_il(1)=0
2070            i1=tdcfl(ch11a(tof11_i),hb11a(tof11_i))
2071            i2=tdcfl(ch11b(tof11_i),hb11b(tof11_i))
2072            if ((i1.eq.1).or.(i2.eq.1)) w_il(1) = 1  ! tdcflag
2073                                                          endif
2074                               endif
2075    
2076            if (tof12_i.gt.0) then
2077            if ((tofmask(ch12a(tof12_i),hb12a(tof12_i)).gt.0).or.
2078         &   (tofmask(ch12b(tof12_i),hb12b(tof12_i)).gt.0)) then
2079            w_il(2)=0
2080            i1=tdcfl(ch12a(tof12_i),hb12a(tof12_i))
2081            i2=tdcfl(ch12b(tof12_i),hb12b(tof12_i))
2082            if ((i1.eq.1).or.(i2.eq.1)) w_il(2) = 1  ! tdcflag
2083                                                          endif
2084                               endif
2085    
2086            if (tof21_i.gt.0) then
2087            if ((tofmask(ch21a(tof21_i),hb21a(tof21_i)).gt.0).or.
2088         &   (tofmask(ch21b(tof21_i),hb21b(tof21_i)).gt.0)) then
2089            w_il(3)=0
2090            i1=tdcfl(ch21a(tof21_i),hb21a(tof21_i))
2091            i2=tdcfl(ch21b(tof21_i),hb21b(tof21_i))
2092            if ((i1.eq.1).or.(i2.eq.1)) w_il(3) = 1  ! tdcflag
2093                                                          endif
2094                               endif
2095    
2096            if (tof22_i.gt.0) then
2097            if ((tofmask(ch22a(tof22_i),hb22a(tof22_i)).gt.0).or.
2098         &   (tofmask(ch22b(tof22_i),hb22b(tof22_i)).gt.0)) then
2099            w_il(4)=0
2100            i1=tdcfl(ch22a(tof22_i),hb22a(tof22_i))
2101            i2=tdcfl(ch22b(tof22_i),hb22b(tof22_i))
2102            if ((i1.eq.1).or.(i2.eq.1)) w_il(4) = 1  ! tdcflag
2103                                                          endif
2104                               endif
2105    
2106            if (tof31_i.gt.0) then
2107            if ((tofmask(ch31a(tof31_i),hb11a(tof31_i)).gt.0).or.
2108         &   (tofmask(ch31b(tof31_i),hb31b(tof31_i)).gt.0)) then
2109            w_il(5)=0
2110            i1=tdcfl(ch31a(tof31_i),hb31a(tof31_i))
2111            i2=tdcfl(ch31b(tof31_i),hb31b(tof31_i))
2112            if ((i1.eq.1).or.(i2.eq.1)) w_il(5) = 1  ! tdcflag
2113                                                          endif
2114                               endif
2115    
2116            if (tof32_i.gt.0) then
2117            if ((tofmask(ch32a(tof32_i),hb32a(tof32_i)).gt.0).or.
2118         &   (tofmask(ch32b(tof32_i),hb32b(tof32_i)).gt.0)) then
2119            w_il(6)=0
2120            i1=tdcfl(ch32a(tof32_i),hb32a(tof32_i))
2121            i2=tdcfl(ch32b(tof32_i),hb32b(tof32_i))
2122            if ((i1.eq.1).or.(i2.eq.1)) w_il(6) = 1  ! tdcflag
2123                                                          endif
2124                               endif
2125    
2126    C------------------------------------------------------------------------
2127    C---  Set weights for the 12 measurements using information for top and bottom:
2128    C---  if no measurements: weight = set to very high value=> not used
2129    C---  top or bottom artificial: weight*sqrt(2)
2130    C---  top and bottom artificial: weight*sqrt(2)*sqrt(2)
2131    
2132           DO jj=1,12
2133           if (jj.le.4)           xhelp = 0.11            ! S1-S3
2134           if ((jj.gt.4).and.(jj.le.8)) xhelp = 0.18      ! S2-S3
2135           if (jj.gt.8)           xhelp = 0.28            ! S1-S2
2136           if ((w_il(itop(jj)).eq.1000.).and.(w_il(ibot(jj)).eq.1000.))
2137         &   xhelp = 1.E09
2138           if ((w_il(itop(jj)).eq.1).or.(w_il(ibot(jj)).eq.1.))
2139         &   xhelp = xhelp*1.414
2140           if ((w_il(itop(jj)).eq.1).and.(w_il(ibot(jj)).eq.1.))
2141         &   xhelp = xhelp*2.
2142           w_i(jj) = 1./xhelp
2143           ENDDO
2144    
2145    C========================================================================
2146    C--- Calculate mean beta for the first time -----------------------------
2147    C--- We are using "1/beta" since its error is gaussian ------------------
2148    
2149            icount=0
2150            sw=0.
2151            sxw=0.
2152            beta_mean=100.
2153    
2154            DO jj=1,12
2155            IF ((abs(1./b(jj)).gt.0.1).and.(abs(1./b(jj)).lt.15.)) THEN
2156                icount = icount+1
2157                sxw    = sxw + (1./b(jj))*w_i(jj)*w_i(jj)
2158                sw     = sw + w_i(jj)*w_i(jj)
2159            ENDIF
2160            ENDDO
2161    
2162            if (icount.gt.0) beta_mean=1./(sxw/sw)
2163            beta_mean_inv = 1./beta_mean
2164    
2165      
2166    C--- Calculate beta for the second time, use residuals of the single
2167    C--- measurements to get a chi2 value
2168    
2169            icount  = 0
2170            sw      = 0.
2171            sxw     = 0.
2172            betachi = 100.
2173            chi2    = 0.
2174            quality = 0.
2175    
2176            DO jj=1,12
2177            IF ((abs(1./b(jj)).gt.0.1).and.(abs(1./b(jj)).lt.15.)
2178         &                                .and.(w_i(jj).GT.0.01)) THEN
2179                res    = beta_mean_inv - (1./b(jj)) ;
2180                if (abs(res*w_i(jj)).lt.resmax) THEN
2181                chi2   = chi2 + (res*w_i(jj))**2.
2182                icount = icount+1
2183                sxw    = sxw + (1./b(jj))*w_i(jj)*w_i(jj)
2184                sw     = sw + w_i(jj)*w_i(jj)
2185                                                ENDIF
2186            ENDIF
2187            ENDDO
2188    
2189    c        quality = sw
2190            quality = sqrt(sw)
2191    
2192            if (icount.eq.0) chi2 = 1000.
2193            if (icount.gt.0) chi2 = chi2/(icount)
2194    
2195            if (icount.gt.0) betachi=1./(sxw/sw);
2196    
2197            beta_mean=100.
2198            if ((chi2.lt.chi2cut).and.(quality.gt.qualitycut))
2199         &  beta_mean = betachi
2200            newbeta = beta_mean
2201    
2202            END
2203    
2204    C****************************************************************************
2205    

Legend:
Removed from v.1.1.1.1  
changed lines
  Added in v.1.13

  ViewVC Help
Powered by ViewVC 1.1.23