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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.13 - (show annotations) (download)
Mon Nov 23 09:50:50 2009 UTC (15 years, 1 month ago) by mocchiut
Branch: MAIN
Changes since 1.12: +5 -182 lines
Cleanup of both C++ and F77 routines

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()
38 c
39 IMPLICIT NONE
40 C
41 include 'input_tof.txt'
42 include 'output_tof.txt'
43 include 'tofcomm.txt'
44
45 INTEGER icounter
46 DATA icounter / 0/
47
48 LOGICAL check
49 REAL secure
50
51 INTEGER j,hitvec(6)
52
53 REAL dx,dy,dr,ds
54 REAL yhelp,yhelp1,yhelp2,xhelp,xhelp1,xhelp2
55 REAL c1,c2
56
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
62 INTEGER tof12_j,tof22_j,tof32_j
63
64 c value for status of each PM-data
65 c first index : 1 = left, 2 = right
66 c second index : 1... number of paddle
67 INTEGER tof11_event(2,8),tof12_event(2,6)
68 INTEGER tof21_event(2,2),tof22_event(2,2)
69 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 theta13
108 C-- DATA ZTOF/53.74,53.04,23.94,23.44,-23.49,-24.34/ !Sergio 9.05.2006
109 REAL tofarm12
110 PARAMETER (tofarm12 = 29.70) ! from 53.39 to 23.69
111 REAL tofarm23
112 PARAMETER (tofarm23 = 47.61) ! from 23.69 to -23.92
113 REAL tofarm13
114 PARAMETER (tofarm13 = 77.31) ! from 53.39 to -23.92
115
116 REAL hepratio
117
118 INTEGER ihelp
119 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---------------------------------------
131 C
132 C Begin !
133 C
134 TOFL2COM = 0
135 C
136 C CALCULATE COMMON VARIABLES
137 C
138 C-------------------------------------------------------------------
139
140 if (ifst.eq.0) then
141
142 ifst=1
143
144 C amplitude has to be 'secure' higher than pedestal for an adc event
145 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
171
172 do i=1,13
173 betatof_a(i) = 100. ! As in "troftrk.for"
174 enddo
175
176 do i=1,6
177 hitvec(i) = -1
178 enddo
179
180 do i=1,4
181 do j=1,12
182 adctof_c(i,j) = 1000.
183 enddo
184 enddo
185
186
187 do i=1,4
188 do j=1,12
189 tdc_c(i,j) = 4095.
190 enddo
191 enddo
192
193
194 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 --------------------------------
226
227 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
231 tof11(1,j,2) = pc_adc(adc(ch11a(j),hb11a(j)))
232 tof11(2,j,2) = pc_adc(adc(ch11b(j),hb11b(j)))
233 tof11(1,j,1) = (tdc(ch11a(j),hb11a(j)))
234 tof11(2,j,1) = (tdc(ch11b(j),hb11b(j)))
235 enddo
236
237
238 do j=1,6
239 tof12(1,j,2) = pc_adc(adc(ch12a(j),hb12a(j)))
240 tof12(2,j,2) = pc_adc(adc(ch12b(j),hb12b(j)))
241 tof12(1,j,1) = (tdc(ch12a(j),hb12a(j)))
242 tof12(2,j,1) = (tdc(ch12b(j),hb12b(j)))
243 enddo
244
245 do j=1,2
246 tof21(1,j,2) = pc_adc(adc(ch21a(j),hb21a(j)))
247 tof21(2,j,2) = pc_adc(adc(ch21b(j),hb21b(j)))
248 tof21(1,j,1) = (tdc(ch21a(j),hb21a(j)))
249 tof21(2,j,1) = (tdc(ch21b(j),hb21b(j)))
250 enddo
251
252 do j=1,2
253 tof22(1,j,2) = pc_adc(adc(ch22a(j),hb22a(j)))
254 tof22(2,j,2) = pc_adc(adc(ch22b(j),hb22b(j)))
255 tof22(1,j,1) = (tdc(ch22a(j),hb22a(j)))
256 tof22(2,j,1) = (tdc(ch22b(j),hb22b(j)))
257 enddo
258
259 do j=1,3
260 tof31(1,j,2) = pc_adc(adc(ch31a(j),hb31a(j)))
261 tof31(2,j,2) = pc_adc(adc(ch31b(j),hb31b(j)))
262 tof31(1,j,1) = (tdc(ch31a(j),hb31a(j)))
263 tof31(2,j,1) = (tdc(ch31b(j),hb31b(j)))
264 enddo
265
266 do j=1,3
267 tof32(1,j,2) = pc_adc(adc(ch32a(j),hb32a(j)))
268 tof32(2,j,2) = pc_adc(adc(ch32b(j),hb32b(j)))
269 tof32(1,j,1) = (tdc(ch32a(j),hb32a(j)))
270 tof32(2,j,1) = (tdc(ch32b(j),hb32b(j)))
271 enddo
272
273 C----------------------------------------------------------------------
274
275 DO i = 1,8
276 if (abs(tof11(1,i,itdc)).gt.10000.) tof11(1,i,itdc)= 10000.
277 if (abs(tof11(2,i,itdc)).gt.10000.) tof11(2,i,itdc)= 10000.
278 if (abs(tof11(1,i,iadc)).gt.10000.) tof11(1,i,iadc)= 10000.
279 if (abs(tof11(2,i,iadc)).gt.10000.) tof11(2,i,iadc)= 10000.
280 ENDDO
281
282 DO i = 1,6
283 if (abs(tof12(1,i,itdc)).gt.10000.) tof12(1,i,itdc)= 10000.
284 if (abs(tof12(2,i,itdc)).gt.10000.) tof12(2,i,itdc)= 10000.
285 if (abs(tof12(1,i,iadc)).gt.10000.) tof12(1,i,iadc)= 10000.
286 if (abs(tof12(2,i,iadc)).gt.10000.) tof12(2,i,iadc)= 10000.
287 ENDDO
288
289
290 DO i = 1,2
291 if (abs(tof21(1,i,itdc)).gt.10000.) tof21(1,i,itdc)= 10000.
292 if (abs(tof21(2,i,itdc)).gt.10000.) tof21(2,i,itdc)= 10000.
293 if (abs(tof21(1,i,iadc)).gt.10000.) tof21(1,i,iadc)= 10000.
294 if (abs(tof21(2,i,iadc)).gt.10000.) tof21(2,i,iadc)= 10000.
295 ENDDO
296
297 DO i = 1,2
298 if (abs(tof22(1,i,itdc)).gt.10000.) tof22(1,i,itdc)= 10000.
299 if (abs(tof22(2,i,itdc)).gt.10000.) tof22(2,i,itdc)= 10000.
300 if (abs(tof22(1,i,iadc)).gt.10000.) tof22(1,i,iadc)= 10000.
301 if (abs(tof22(2,i,iadc)).gt.10000.) tof22(2,i,iadc)= 10000.
302 ENDDO
303
304 DO i = 1,3
305 if (abs(tof31(1,i,itdc)).gt.10000.) tof31(1,i,itdc)= 10000.
306 if (abs(tof31(2,i,itdc)).gt.10000.) tof31(2,i,itdc)= 10000.
307 if (abs(tof31(1,i,iadc)).gt.10000.) tof31(1,i,iadc)= 10000.
308 if (abs(tof31(2,i,iadc)).gt.10000.) tof31(2,i,iadc)= 10000.
309 ENDDO
310
311 DO i = 1,3
312 if (abs(tof32(1,i,itdc)).gt.10000.) tof32(1,i,itdc)= 10000.
313 if (abs(tof32(2,i,itdc)).gt.10000.) tof32(2,i,itdc)= 10000.
314 if (abs(tof32(1,i,iadc)).gt.10000.) tof32(1,i,iadc)= 10000.
315 if (abs(tof32(2,i,iadc)).gt.10000.) tof32(2,i,iadc)= 10000.
316 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----------------------------------------------------------------
360 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----------------------------------------------------------------
379
380 C upper tof S11
381 DO i = 1,8
382
383 DO j = 1,2
384 tof11_event(j,i) = none_ev
385 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
387 ENDDO
388 ENDDO
389
390 c find single paddle in upper tof with tdc and adc signal
391 tof11_i = none_find
392 tof11_j = none_find
393 check = .TRUE.
394 DO i = 1, 8
395 IF ((tof11_event(left,i).GE.1).AND.(tof11_event(right,i).GE.1))
396 + THEN
397 c check if an other paddle has also an event - then set flag
398 tof11_j = tof11_j + 2**(i-1)
399 IF (check.EQV..TRUE.) THEN
400 IF (tof11_i.EQ.none_find) THEN
401 tof11_i = i
402 ELSE
403 tof11_i = -1
404 check = .FALSE.
405 ENDIF
406 ENDIF
407 ENDIF
408 ENDDO
409
410
411 C upper tof S12
412 DO i = 1,6
413 DO j = 1,2
414 tof12_event(j,i) = none_ev
415 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
417 ENDDO
418 ENDDO
419
420 c find single paddle in upper tof with tdc and adc signal
421 tof12_i = none_find
422 tof12_j = none_find
423 check = .TRUE.
424 DO i = 1, 6
425 IF ((tof12_event(left,i).GE.1).AND.(tof12_event(right,i).GE.1))
426 + THEN
427 c check if an other paddle has also an event - then set flag
428 tof12_j = tof12_j + 2**(i-1)
429 IF (check.EQV..TRUE.) THEN
430 IF (tof12_i.EQ.none_find) THEN
431 tof12_i = i
432 ELSE
433 tof12_i = -1
434 check = .FALSE.
435 ENDIF
436 ENDIF
437 ENDIF
438 ENDDO
439
440
441 C middle tof S21
442 DO i = 1,2
443 DO j = 1,2
444 tof21_event(j,i) = none_ev
445 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
447 ENDDO
448 ENDDO
449
450 c find single paddle in upper tof with tdc and adc signal
451 tof21_i = none_find
452 tof21_j = none_find
453 check = .TRUE.
454 DO i = 1, 2
455 IF ((tof21_event(left,i).GE.1).AND.(tof21_event(right,i).GE.1))
456 + THEN
457 c check if an other paddle has also an event - then set flag
458 tof21_j = tof21_j + 2**(i-1)
459 IF (check.EQV..TRUE.) THEN
460 IF (tof21_i.EQ.none_find) THEN
461 tof21_i = i
462 ELSE
463 tof21_i = -1
464 check = .FALSE.
465 ENDIF
466 ENDIF
467 ENDIF
468 ENDDO
469
470 C middle tof S22
471 DO i = 1,2
472 DO j = 1,2
473 tof22_event(j,i) = none_ev
474 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
476 ENDDO
477 ENDDO
478
479 c find single paddle in upper tof with tdc and adc signal
480 tof22_i = none_find
481 tof22_j = none_find
482 check = .TRUE.
483 DO i = 1, 2
484 IF ((tof22_event(left,i).GE.1).AND.(tof22_event(right,i).GE.1))
485 + THEN
486 c check if an other paddle has also an event - then set flag
487 tof22_j = tof22_j + 2**(i-1)
488 IF (check.EQV..TRUE.) THEN
489 IF (tof22_i.EQ.none_find) THEN
490 tof22_i = i
491 ELSE
492 tof22_i = -1
493 check = .FALSE.
494 ENDIF
495 ENDIF
496 ENDIF
497 ENDDO
498
499
500 C bottom tof S31
501 DO i = 1,3
502 DO j = 1,2
503 tof31_event(j,i) = none_ev
504 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
506 ENDDO
507 ENDDO
508
509 c find single paddle in upper tof with tdc and adc signal
510 tof31_i = none_find
511 tof31_j = none_find
512 check = .TRUE.
513 DO i = 1, 3
514 IF ((tof31_event(left,i).GE.1).AND.(tof31_event(right,i).GE.1))
515 + THEN
516 c check if an other paddle has also an event - then set flag
517 tof31_j = tof31_j + 2**(i-1)
518 IF (check.EQV..TRUE.) THEN
519 IF (tof31_i.EQ.none_find) THEN
520 tof31_i = i
521 ELSE
522 tof31_i = -1
523 check = .FALSE.
524 ENDIF
525 ENDIF
526 ENDIF
527 ENDDO
528
529 C bottom tof S32
530 DO i = 1,3
531 DO j = 1,2
532 tof32_event(j,i) = none_ev
533 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
535 ENDDO
536 ENDDO
537
538 c find single paddle in upper tof with tdc and adc signal
539 tof32_i = none_find
540 tof32_j = none_find
541 check = .TRUE.
542 DO i = 1, 3
543 IF ((tof32_event(left,i).GE.1).AND.(tof32_event(right,i).GE.1))
544 + THEN
545 c check if an other paddle has also an event - then set flag
546 tof32_j = tof32_j + 2**(i-1)
547 IF (check.EQV..TRUE.) THEN
548 IF (tof32_i.EQ.none_find) THEN
549 tof32_i = i
550 ELSE
551 tof32_i = -1
552 check = .FALSE.
553 ENDIF
554 ENDIF
555 ENDIF
556 ENDDO
557
558 do i=1,6
559 tof_i_flag(i)=0
560 tof_j_flag(i)=0
561 enddo
562
563 tof_i_flag(1)=tof11_i
564 tof_i_flag(2)=tof12_i
565 tof_i_flag(3)=tof21_i
566 tof_i_flag(4)=tof22_i
567 tof_i_flag(5)=tof31_i
568 tof_i_flag(6)=tof32_i
569
570 tof_j_flag(1)=tof11_j
571 tof_j_flag(2)=tof12_j
572 tof_j_flag(3)=tof21_j
573 tof_j_flag(4)=tof22_j
574 tof_j_flag(5)=tof31_j
575 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
585 C------------------------------------------------------------------
586 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------------------------------------------------------------------
592
593 do i=1,3
594 xtofpos(i)=100.
595 ytofpos(i)=100.
596 enddo
597
598 C-----------------------------S1 --------------------------------
599
600 IF (tof11_i.GT.none_find) THEN
601 ytofpos(1) = ((tof11(1,tof11_i,itdc)-tof11(2,tof11_i,itdc))/2.
602 + -y_coor_lin11c(tof11_i,offset))/y_coor_lin11c(tof11_i,slope)
603 endif
604
605 IF (tof12_i.GT.none_find) THEN
606 xtofpos(1) = ((tof12(1,tof12_i,itdc)-tof12(2,tof12_i,itdc))/2.
607 + -x_coor_lin12c(tof12_i,offset))/x_coor_lin12c(tof12_i,slope)
608 endif
609
610
611 C-----------------------------S2 --------------------------------
612
613 IF (tof21_i.GT.none_find) THEN
614 xtofpos(2) = ((tof21(1,tof21_i,itdc)-tof21(2,tof21_i,itdc))/2.
615 + -x_coor_lin21c(tof21_i,offset))/x_coor_lin21c(tof21_i,slope)
616 endif
617
618 IF (tof22_i.GT.none_find) THEN
619 ytofpos(2) = ((tof22(1,tof22_i,itdc)-tof22(2,tof22_i,itdc))/2.
620 + -y_coor_lin22c(tof22_i,offset))/y_coor_lin22c(tof22_i,slope)
621 endif
622
623
624 C-----------------------------S3 --------------------------------
625
626 IF (tof31_i.GT.none_find) THEN
627 ytofpos(3) = ((tof31(1,tof31_i,itdc)-tof31(2,tof31_i,itdc))/2.
628 + -y_coor_lin31c(tof31_i,offset))/y_coor_lin31c(tof31_i,slope)
629 endif
630
631 IF (tof32_i.GT.none_find) THEN
632 xtofpos(3) = ((tof32(1,tof32_i,itdc)-tof32(2,tof32_i,itdc))/2.
633 + -x_coor_lin32c(tof32_i,offset))/x_coor_lin32c(tof32_i,slope)
634 endif
635
636
637 C----------------------------------------------------------------------
638 C--------------------- zenith angle theta ---------------------------
639 C----------------------------------------------------------------------
640
641 xhelp1=0.
642 if (tof11_i.GT.none_find) xhelp1=tof11_x(tof11_i)
643 if (xtofpos(1).lt.100) xhelp1=xtofpos(1)
644
645 yhelp1=0.
646 if (tof12_i.GT.none_find) yhelp1=tof12_y(tof12_i)
647 if (ytofpos(1).lt.100) yhelp1=ytofpos(1)
648
649
650 yhelp2=0.
651 if (tof32_i.GT.none_find) yhelp2=tof32_y(tof32_i)
652 if (ytofpos(3).lt.100) yhelp2=ytofpos(3)
653
654 xhelp2=0.
655 if (tof31_i.GT.none_find) xhelp2=tof31_x(tof31_i)
656 if (xtofpos(3).lt.100) xhelp2=xtofpos(3)
657
658
659 dx=0.
660 dy=0.
661 dr=0.
662 theta13 = 0.
663
664 dx = xhelp1 - xhelp2
665 dy = yhelp1 - yhelp2
666 dr = sqrt(dx*dx+dy*dy)
667 theta13 = atan(dr/tofarm13)
668
669
670 C----------------------------------------------------------------------
671 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----------------------------------------------------------------------
677
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
686 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/
688 c DATA tof21_y/ 3.75,-3.75/ ! paddles in different order
689 c DATA tof22_x/ -4.5,4.5/
690 c DATA tof31_x/ -6.0,0.,6.0/
691 c DATA tof32_y/ -5.0,0.0,5.0/
692
693
694 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
703 if (adc(ch11a(i),hb11a(i)).eq.4095) then
704 xkorr = atten(left,11,i,yhelp)
705 if (iz.le.1) xkorr=xkorr/hepratio
706 tof11(left,i,iadc)=xkorr/cos(theta13)
707 adcflagtof(ch11a(i),hb11a(i)) = 1
708 endif
709 if (adc(ch11b(i),hb11b(i)).eq.4095) then
710 xkorr = atten(right,11,i,yhelp)
711 if (iz.le.1) xkorr=xkorr/hepratio
712 tof11(right,i,iadc)=xkorr/cos(theta13)
713 adcflagtof(ch11b(i),hb11b(i)) = 1
714 endif
715 ENDIF
716
717 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
724 if (adc(ch12a(i),hb12a(i)).eq.4095) then
725 xkorr = atten(left,12,i,xhelp)
726 if (iz.le.1) xkorr=xkorr/hepratio
727 tof12(left,i,iadc) = xkorr/cos(theta13)
728 adcflagtof(ch12a(i),hb12a(i)) = 1
729 endif
730 if (adc(ch12b(i),hb12b(i)).eq.4095) then
731 xkorr = atten(right,12,i,xhelp)
732 if (iz.le.1) xkorr=xkorr/hepratio
733 tof12(right,i,iadc) = xkorr/cos(theta13)
734 adcflagtof(ch12b(i),hb12b(i)) = 1
735 endif
736 ENDIF
737
738 C-----------------------------S2 --------------------------------
739
740 c xhelp=0.
741 xhelp=100. ! WM
742 if (tof22_i.GT.none_find) xhelp=tof22_x(tof22_i)
743 if (xtofpos(2).lt.100) xhelp=xtofpos(2)
744
745 IF (tof21_i.GT.none_find.AND.abs(xhelp).lt.100) THEN
746 i = tof21_i
747 if (adc(ch21a(i),hb21a(i)).eq.4095) then
748 xkorr = atten(left,21,i,xhelp)
749 if (iz.le.1) xkorr=xkorr/hepratio
750 tof21(left,i,iadc) = xkorr/cos(theta13)
751 adcflagtof(ch21a(i),hb21a(i)) = 1
752 endif
753 if (adc(ch21b(i),hb21b(i)).eq.4095) then
754 xkorr = atten(right,21,i,xhelp)
755 if (iz.le.1) xkorr=xkorr/hepratio
756 tof21(right,i,iadc) = xkorr/cos(theta13)
757 adcflagtof(ch21b(i),hb21b(i)) = 1
758 endif
759 ENDIF
760
761
762 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
769 if (adc(ch22a(i),hb22a(i)).eq.4095) then
770 xkorr = atten(left,22,i,yhelp)
771 if (iz.le.1) xkorr=xkorr/hepratio
772 tof22(left,i,iadc) = xkorr/cos(theta13)
773 adcflagtof(ch22a(i),hb22a(i)) = 1
774 endif
775 if (adc(ch22b(i),hb22b(i)).eq.4095) then
776 xkorr = atten(right,22,i,yhelp)
777 if (iz.le.1) xkorr=xkorr/hepratio
778 tof22(right,i,iadc) = xkorr/cos(theta13)
779 adcflagtof(ch22b(i),hb22b(i)) = 1
780 endif
781 ENDIF
782
783 C-----------------------------S3 --------------------------------
784
785 c yhelp=0.
786 yhelp=100. ! WM
787 if (tof32_i.GT.none_find) yhelp=tof32_y(tof32_i)
788 if (ytofpos(3).lt.100) yhelp=ytofpos(3)
789
790 IF (tof31_i.GT.none_find.AND.abs(yhelp).lt.100) THEN
791 i = tof31_i
792 if (adc(ch31a(i),hb31a(i)).eq.4095) then
793 xkorr = atten(left,31,i,yhelp)
794 if (iz.le.1) xkorr=xkorr/hepratio
795 tof31(left,i,iadc) = xkorr/cos(theta13)
796 adcflagtof(ch31a(i),hb31a(i)) = 1
797 endif
798 if (adc(ch31b(i),hb31b(i)).eq.4095) then
799 xkorr = atten(right,31,i,yhelp)
800 if (iz.le.1) xkorr=xkorr/hepratio
801 tof31(right,i,iadc) = xkorr/cos(theta13)
802 adcflagtof(ch31b(i),hb31b(i)) = 1
803 endif
804 ENDIF
805
806 c xhelp=0.
807 xhelp=100. ! WM
808 if (tof31_i.GT.none_find) xhelp=tof31_x(tof31_i)
809 if (xtofpos(3).lt.100) xhelp=xtofpos(3)
810
811 IF (tof32_i.GT.none_find.AND.abs(xhelp).lt.100) THEN
812 i = tof32_i
813 if (adc(ch32a(i),hb32a(i)).eq.4095) then
814 xkorr = atten(left,32,i,xhelp)
815 if (iz.le.1) xkorr=xkorr/hepratio
816 tof32(left,i,iadc) = xkorr/cos(theta13)
817 adcflagtof(ch32a(i),hb32a(i)) = 1
818 endif
819 if (adc(ch32b(i),hb32b(i)).eq.4095) then
820 xkorr = atten(right,32,i,xhelp)
821 if (iz.le.1) xkorr=xkorr/hepratio
822 tof32(right,i,iadc) = xkorr/cos(theta13)
823 adcflagtof(ch32b(i),hb32b(i)) = 1
824 endif
825 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
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 ------------------------------
1066 C--------------------------------------------------------------------
1067 C-------------------difference of sums -----------------------------
1068 C
1069 C DS = (t1+t2) - t3+t4)
1070 C DS = c1 + c2/beta*cos(theta)
1071 C c2 = 2d/c gives c2 = 2d/(c*TDCresolution) TDC=50ps/channel
1072 C => c2 = ca.60 for 0.45 m c2 = ca.109 for 0.81 m
1073 C since TDC resolution varies slightly c2 has to be calibrated
1074
1075 C S11 - S31
1076
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)
1080 xhelp2 = tof31(1,tof31_i,itdc)+tof31(2,tof31_i,itdc)
1081 ds = xhelp1-xhelp2
1082 ihelp=(tof11_i-1)*3+tof31_i
1083 c1 = k_S11S31(1,ihelp)
1084 if (iz.gt.2) c1 = c1 + k1corrA1
1085 c2 = k_S11S31(2,ihelp)
1086 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
1103
1104 C S11 - S32
1105
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)
1109 xhelp2 = tof32(1,tof32_i,itdc)+tof32(2,tof32_i,itdc)
1110 ds = xhelp1-xhelp2
1111 ihelp=(tof11_i-1)*3+tof32_i
1112 c1 = k_S11S32(1,ihelp)
1113 if (iz.gt.2) c1 = c1 + k1corrA1
1114 c2 = k_S11S32(2,ihelp)
1115 betatof_a(2) = c2/(cos(theta13)*(ds-c1))
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
1134
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)
1138 xhelp2 = tof31(1,tof31_i,itdc)+tof31(2,tof31_i,itdc)
1139 ds = xhelp1-xhelp2
1140 ihelp=(tof12_i-1)*3+tof31_i
1141 c1 = k_S12S31(1,ihelp)
1142 if (iz.gt.2) c1 = c1 + k1corrA1
1143 c2 = k_S12S31(2,ihelp)
1144 betatof_a(3) = c2/(cos(theta13)*(ds-c1))
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
1163
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)
1167 xhelp2 = tof32(1,tof32_i,itdc)+tof32(2,tof32_i,itdc)
1168 ds = xhelp1-xhelp2
1169 ihelp=(tof12_i-1)*3+tof32_i
1170 c1 = k_S12S32(1,ihelp)
1171 if (iz.gt.2) c1 = c1 + k1corrA1
1172 c2 = k_S12S32(2,ihelp)
1173 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
1190
1191 C S21 - S31
1192
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)
1196 xhelp2 = tof31(1,tof31_i,itdc)+tof31(2,tof31_i,itdc)
1197 ds = xhelp1-xhelp2
1198 ihelp=(tof21_i-1)*3+tof31_i
1199 c1 = k_S21S31(1,ihelp)
1200 if (iz.gt.2) c1 = c1 + k1corrB1
1201 c2 = k_S21S31(2,ihelp)
1202 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
1219
1220 C S21 - S32
1221
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)
1225 xhelp2 = tof32(1,tof32_i,itdc)+tof32(2,tof32_i,itdc)
1226 ds = xhelp1-xhelp2
1227 ihelp=(tof21_i-1)*3+tof32_i
1228 c1 = k_S21S32(1,ihelp)
1229 if (iz.gt.2) c1 = c1 + k1corrB1
1230 c2 = k_S21S32(2,ihelp)
1231 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
1248
1249 C S22 - S31
1250
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)
1254 xhelp2 = tof31(1,tof31_i,itdc)+tof31(2,tof31_i,itdc)
1255 ds = xhelp1-xhelp2
1256 ihelp=(tof22_i-1)*3+tof31_i
1257 c1 = k_S22S31(1,ihelp)
1258 if (iz.gt.2) c1 = c1 + k1corrB1
1259 c2 = k_S22S31(2,ihelp)
1260 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
1277
1278 C S22 - S32
1279
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)
1283 xhelp2 = tof32(1,tof32_i,itdc)+tof32(2,tof32_i,itdc)
1284 ds = xhelp1-xhelp2
1285 ihelp=(tof22_i-1)*3+tof32_i
1286 c1 = k_S22S32(1,ihelp)
1287 if (iz.gt.2) c1 = c1 + k1corrB1
1288 c2 = k_S22S32(2,ihelp)
1289 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
1306
1307 C S11 - S21
1308
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)
1312 xhelp2 = tof21(1,tof21_i,itdc)+tof21(2,tof21_i,itdc)
1313 ds = xhelp1-xhelp2
1314 ihelp=(tof11_i-1)*2+tof21_i
1315 c1 = k_S11S21(1,ihelp)
1316 if (iz.gt.2) c1 = c1 + k1corrC1
1317 c2 = k_S11S21(2,ihelp)
1318 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
1335
1336 C S11 - S22
1337
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)
1341 xhelp2 = tof22(1,tof22_i,itdc)+tof22(2,tof22_i,itdc)
1342 ds = xhelp1-xhelp2
1343 ihelp=(tof11_i-1)*2+tof22_i
1344 c1 = k_S11S22(1,ihelp)
1345 if (iz.gt.2) c1 = c1 + k1corrC1
1346 c2 = k_S11S22(2,ihelp)
1347 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
1364
1365 C S12 - S21
1366
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)
1370 xhelp2 = tof21(1,tof21_i,itdc)+tof21(2,tof21_i,itdc)
1371 ds = xhelp1-xhelp2
1372 ihelp=(tof12_i-1)*2+tof21_i
1373 c1 = k_S12S21(1,ihelp)
1374 if (iz.gt.2) c1 = c1 + k1corrC1
1375 c2 = k_S12S21(2,ihelp)
1376 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
1393
1394 C S12 - S22
1395
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)
1399 xhelp2 = tof22(1,tof22_i,itdc)+tof22(2,tof22_i,itdc)
1400 ds = xhelp1-xhelp2
1401 ihelp=(tof12_i-1)*2+tof22_i
1402 c1 = k_S12S22(1,ihelp)
1403 if (iz.gt.2) c1 = c1 + k1corrC1
1404 c2 = k_S12S22(2,ihelp)
1405 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
1422
1423 C---------------------------------------------------------
1424 C
1425 C icount=0
1426 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
1448 btemp(i) = betatof_a(i)
1449 enddo
1450
1451 betatof_a(13)=newbeta(1,btemp,hitvec,10.,10.,20.)
1452
1453 C--------------------------------------------------------------
1454 C write(*,*) betatof_a
1455 c write(*,*) xtofpos
1456 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 c 100 continue
1466 continue
1467
1468 C
1469 RETURN
1470 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

  ViewVC Help
Powered by ViewVC 1.1.23