MACRO init path=/shift/shift19/data-d2/a/attal

chain -res
nt/chain res [path]/sync_cscdata_0103.ntuple _
[path]/sync_cscdata_0104.ntuple
 

cd //res
uwfunc //res/10 looping.inc

cd //res

RETURN
 
 

MACRO find_tracks

nt/create 20 'Alons ntuple' 15 ! ! _
evnt_num _
h_slope _
h_intrcp _
v_slope _
v_intrcp _
h_proj _
v_proj _
h_strip _
an_wire _
clct_hs _
alct_wir _
hs_diff _
an_diff _
pt _
pattern
 

exe init

Application Comis Quit
Real Function dump_tele
implicit none
include 'looping.inc'

real
. sum_E,
. sum_Q,
. sum_N,
. sum_D

integer
. event

integer
. max_layer_hits,
. n_layers,
. max_peaks,
. minus,
. plus,
. max_allowed_tracks,
. vertical_strips,
. horizontal_strips
parameter
.(max_layer_hits = 50,
. n_layers = 4,
. max_peaks = 20,
. minus = 1,
. plus = 2,
. max_alowed_tracks = 16,
. vertical_strips = 1,
. horizontal_strips = 2)

integer
. L,H,N,P,D,
. check_sum,
. P1,P2,P3,P4

integer
. hits( 4),
. signal(50,4),
. noise( 50,4),
. strip( 50,4)
 

integer
. peaks( 4), ! # of peaks in layer
. peak_strip( 20,4), ! used as an index to the hit
. neighbors(1:2,20,4)

real
. CoG(20,4),
. e_CoG(20,4)

integer
. n_trk_4(1:2)
real
. trk( 4,
. 16,
. 1:2),
. e_trk( 4,
. 16,
. 1:2)
integer
. trk_peak(4,
. 16,
. 1:2),
. trk_layer(4,
. 16,
. 1:2)

real
. x(4,1:2),
. par_v(2),
. e_par_v(2),
. chi2_v,
. par_h(2),
. e_par_h(2),
. chi2_h,
. dummy

real
. delta_12,
. delta_23,
. delta_34

real
. max_two_layer_diff

parameter
.(max_two_layer_diff = 30.)

real
. delta_delta,
. delta_delta_check,
. delta_12_check,
. delta_23_check,
. delta_34_check,
. trk_c(4)
integer
. P1c,P2c,P3c,P4c

logical
. copied_hits
integer
. copied_hit_candidates
logical
. copied(20,4)

integer
. L1,L2,L3

integer
. T1,T2

real
. x_V(4),
. x_H(4),
. x3_V(3),
. x3_H(3)

real
. tk(4),
. etk(4),
. Slope_V,
. Intercept_V,
. Error_slope_V,
. Error_intercept_V,
. Linear_correlation_constant_V,
. Least_squares_by_DOF_V,
. Slope_H,
. Intercept_H,
. Error_slope_H,
. Error_intercept_H,
. Linear_correlation_constant_H,
. Least_squares_by_DOF_H
 

real
. Projection_H,
. Projection_V,
. Half_strip,
. Anode_wire,
. Err_Proj_H,
. HS_diff,
. Anode_diff,
. trans_mom

real
. Travel_dist
parameter
. (Travel_dist = 50)
* Travel_dist is 2402(dist from center of chamber to center of telescope in cm)

real
. Half_strip_width
parameter
. (Half_strip_width = 0.4917)
* 1/2 strip width at particular beam spot in cm

real
. Wire_spacing
parameter
. (Wire_spacing = 2.54)
* Distance between two anode wires in cm

real
. C_LCT_offset
parameter
. (C_LCT_offset = -1.55)
* dist between zero of silicon beam tele to zero of CLCT board in cm

real
. A_LCT_offset
parameter
. (A_LCT_offset = 0)

real
. Slope_coefficient
parameter
. (Slope_coefficient = 0.1)
 

real
. temp,
. correction,
. temp_offset,
. strip_offset,
. h_dist

integer
. pattern

integer
. n_trk_3(1:2)

real
. trk3( 3,
. 16,
. 1:2),
. e_trk3( 3,
. 16,
. 1:2)
integer
. trk_peak3( 3,
. 16,
. 1:2),
. trk_layer3(3,
. 16,
. 1:2)

logical
. found_layer(4)

real
. ntuple(8)

real
. final_t(8),
. e_final_t(8)

real
. x_zero,
. x_P2
parameter
.(x_zero = 526.,
. x_P2 = 1770.)

real
. offset(4,1:2)

real
. tilt(4,1:2)

offset(1,2) = -6.094
offset(2,2) = 14.868
offset(3,2) = -7.963
offset(4,2) = 0.424

offset(1,1) = 4.528
offset(2,1) = -6.050
offset(3,1) = -1.970
offset(4,1) = 3.506
 

tilt(1,1) = 0.
tilt(2,1) = 0.
tilt(3,1) = -1.5
tilt(4,1) = 0.

tilt(1,2) = -10.5
tilt(2,2) = 2.25
tilt(3,2) = 3.75
tilt(4,2) = 4.4
 

x(1,1) = 245. - 526.
x(2,1) = 427. - 526.
x(3,1) = 587. - 526.
x(4,1) = 754. - 526.

x(1,2) = 283. - 526.
x(2,2) = 466. - 526.
x(3,2) = 626. - 526.
x(4,2) = 792. - 526.

x_V(1) = 245. - 526.
x_V(2) = 427. - 526.
x_V(3) = 587. - 526.
x_V(4) = 754. - 526.

x_H(1) = 283. - 526.
x_H(2) = 466. - 526.
x_H(3) = 626. - 526.
x_H(4) = 792. - 526.
 

event = event_number
 

********************* Selects tracks ************************
 

do D = 1,2

do L = 1,4
hits(L) = 0
enddo

do H = 1,n_si_tele

if(mod(si_tele_layer(H),2) .eq. 0 .and.
D .eq. 2) then
L = si_tele_layer(H)/2
hits(L) = hits(L) + 1

if( hits(L) .gt. 50 ) then
* write(*,*) 'ERROR: Increase # of max layer hits',hits(L)
goto 999
endif

strip(hits(L),L) = si_tele_strip(H)
noise(hits(L),L) = si_tele_noise(H)

if( si_tele_signal(H) .gt. 32768) then
signal(hits(L),L) = si_tele_signal(H) - 65536
else
signal(hits(L),L) = si_tele_signal(H)
endif
elseif(mod(si_tele_layer(H),2) .eq. 1 .and.
D .eq. 1) then
L = (si_tele_layer(H)+1)/2
hits(L) = hits(L) + 1

if( hits(L) .gt. 50 ) then
* write(*,*) 'ERROR: Increase # of max layer hits',hits(L)
goto 999
endif

strip(hits(L),L) = si_tele_strip(H)
noise(hits(L),L) = si_tele_noise(H)

if( si_tele_signal(H) .gt. 32768) then
signal(hits(L),L) = si_tele_signal(H) - 65536
else
signal(hits(L),L) = si_tele_signal(H)
endif

endif
 

enddo
 

do L = 1,4
do H = 1,hits(L)-1
if(strip(H,L).ge.strip(H+1,L) ) then
* write(*,*) 'ERROR, ordering problem, skipping event.'
* write(*,*) ' Layer, strip n,strip n+1:',L,strip(H, L),
* . strip(H+1,L)
goto 999
endif
enddo
enddo

* Make all of the clusters in a layer:

do L = 1,4
peaks(L) = 0
do H = 1,20
do N = 1,2
neighbors(N,H,L) = 0
enddo
enddo
enddo

do L = 1,4

if( hits(L) .eq. 1) then
peaks(L) = 1
peak_strip( peaks(L),L) = 1
neighbors(1,peaks(L),L) = 0
neighbors( 2,peaks(L),L) = 0
 

elseif(hits(L) .eq. 2) then

if( strip(1,L) .eq. strip(2,L)-1) then
if(signal(1,L) .ge. signal(2,L) ) then

peaks(L) = 1
peak_strip( peaks(L),L) = 1
neighbors(1,peaks(L),L) = 0
neighbors( 2,peaks(L),L) = 1
else
peaks(L) = 1
peak_strip( peaks(L),L) = 2
neighbors(1,peaks(L),L) = 1
neighbors( 2,peaks(L),L) = 0
endif
else
peaks(L) = 1
peak_strip( peaks(L),L) = 1
neighbors(1,peaks(L),L) = 0
neighbors( 2,peaks(L),L) = 0
peaks(L) = 2
peak_strip( peaks(L),L) = 2
neighbors(1,peaks(L),L) = 0
neighbors( 2,peaks(L),L) = 0
endif
 

elseif(hits(L) .ge. 3) then ! first and last strip treated special

if( strip(1,L) .eq. strip(2,L)-1 ) then ! first strip
if( signal(1,L) .ge. signal(2,L) ) then

peaks(L) = peaks(L) + 1
peak_strip( peaks(L),L) = 1
neighbors(1,peaks(L),L) = 0

do N = 1,hits(L)-1
if(signal(N,L) .ge. signal(N+1,L) .and.
. strip(N,L) .eq. strip(N+1,L)-1 ) then

neighbors( 2,peaks(L),L) = neighbors( 2,peaks(L),L) +1

else
goto 11
endif
enddo
11 continue
endif
else
peaks(L) = peaks(L) + 1
peak_strip(peaks(L),L) = 1
neighbors(1,peaks(L),L) = 0
neighbors( 2,peaks(L),L) = 0
endif
 
 

do H = 2,hits(L)-1

if( strip(H,L) .eq. strip(H-1,L)+1 .and.
. strip(H,L) .eq. strip(H+1,L)-1 ) then
if( signal(H,L) .gt. signal(H-1,L) .and.
. signal(H,L) .ge. signal(H+1,L) ) then

peaks(L) = peaks(L) + 1
peak_strip( peaks(L),L) = H

do N = H,hits(L)-1
if(signal(N,L) .ge. signal(N+1,L) .and.
. strip(N,L) .eq. strip(N+1,L)-1 ) then

neighbors( 2,peaks(L),L) =neighbors( 2,peaks(L),L)+1

else
goto 13
endif
enddo
13 continue

do N = H,2,-1
if(signal(N,L) .gt. signal(N-1,L) .and.
. strip(N,L) .eq. strip(N-1,L)+1 ) then

neighbors(1,peaks(L),L) =neighbors(1,peaks(L),L)+1

else
goto 14
endif
enddo
14 continue

endif

elseif(strip( H,L) .eq. strip( H-1,L)+1 .and.
. strip( H,L) .ne. strip( H+1,L)-1 ) then
if( signal(H,L) .gt. signal(H-1,L) ) then

peaks(L) = peaks(L) + 1
peak_strip( peaks(L),L) = H
neighbors( 2,peaks(L),L) = 0

do N = H,2,-1
if(signal(N,L) .gt. signal(N-1,L) .and.
. strip(N,L) .eq. strip(N-1,L)+1 ) then

neighbors(1,peaks(L),L) =neighbors(1,peaks(L),L)+1

else
goto 15
endif
enddo
15 continue
endif

elseif(strip( H,L) .ne. strip( H-1,L)+1 .and.
. strip( H,L) .eq. strip( H+1,L)-1 ) then
if( signal(H,L) .ge. signal(H+1,L) ) then

peaks(L) = peaks(L) + 1
peak_strip( peaks(L),L) = H
neighbors(1,peaks(L),L) = 0

do N = H,hits(L)-1
if(signal(N,L) .ge. signal(N+1,L) .and.
. strip(N,L) .eq. strip(N+1,L)-1 ) then

neighbors( 2,peaks(L),L) =neighbors( 2,peaks(L),L)+1

else
goto 16
endif
enddo
16 continue

endif

elseif(strip( H,L) .ne. strip( H-1,L)+1 .and.
. strip( H,L) .ne. strip( H+1,L)-1 ) then

peaks(L) = peaks(L) + 1
peak_strip( peaks(L),L) = H
neighbors(1,peaks(L),L) = 0
neighbors( 2,peaks(L),L) = 0

endif
if(peaks(L) .ge. 20) then
* write(*,*)'ERROR: increase number of max peaks!',peaks(L)
goto 999
endif
enddo
 

if( strip(hits(L),L) .eq. strip(hits(L)-1,L)+1) then ! last strip
if( signal(hits(L),L) .gt. signal(hits(L)-1,L)) then

peaks(L) = peaks(L) + 1
peak_strip( peaks(L),L) = hits(L)
neighbors(2,peaks(L),L) = 0

do N = hits(L),2,-1
if(signal(N,L) .gt. signal(N-1,L) .and.
. strip(N,L) .eq. strip(N-1,L)+1 ) then

neighbors(1,peaks(L),L) = neighbors(1,peaks(L),L) +1

else
goto 12
endif
enddo
12 continue
endif
else
peaks(L) = peaks(L) + 1
peak_strip( peaks(L),L) = hits(L)
neighbors(1,peaks(L),L) = 0
neighbors( 2,peaks(L),L) = 0
endif

endif
enddo

* OK, now we resolve the conflicts by splitting charge or by assigning to closest.

do L = 1,4
do P = 1,peaks(L)-1

if( strip(peak_strip( P,L),L) .ge.
. strip(peak_strip(P+1,L),L) ) then
write(*,*) 'ERROR: strip ordering wrong',
. strip(peak_strip( P,L),L),strip(peak_strip(P+1,L),L)
goto 999
endif
 

if( strip(peak_strip( P,L),L) + neighbors( 2, P,L) .ge.
. strip(peak_strip(P+1,L),L) - neighbors(1,P+1,L) ) then

if(1) then
* write(*,*) 'WARNING, found conflict :',
* . strip(peak_strip( P,L),L) + neighbors( 2, P,L),
* . strip(peak_strip(P+1,L),L) - neighbors(1,P+1,L)
endif
 

if( neighbors( 2, P,L) .eq.
. neighbors(1,P+1,L) ) then

if(1) then
* write(*,*)'SOLUTION: split the charge in half:',
* . neighbors( 2, P,L)
endif

do N = P+1,peaks(L)
peak_strip(N,L) = peak_strip(N,L) + 1
enddo

hits(L) = hits(L) + 1

if(hits(L) .gt. max_layer_hits ) then
* write(*,*) 'ERROR: can''t add more hits: ', hits(L)
goto 999
endif

do H = hits(L)-1,peak_strip(P,L)+neighbors( 2, P,L)+1,-1
signal(H+1,L) = signal(H,L)
noise( H+1,L) = noise( H,L)
strip( H+1,L) = strip( H,L)
enddo

H = peak_strip(P,L)+neighbors( 2, P,L)

noise( H+1,L) = noise( H,L)
strip( H+1,L) = strip( H,L)
 

signal(H+1,L) = float(
.
. signal(H,L)*
. signal(peak_strip(P+1,L), L)/
. (signal(peak_strip( P,L), L)+signal(peak_strip(P+1,L),L)))

signal(H,L) = float(
.
. signal(H,L)*
. signal(peak_strip( P,L), L)/
. (signal(peak_strip( P,L), L)+signal(peak_strip(P+1,L),L)))
 

elseif( neighbors( 2, P,L) .gt.
. neighbors(1,P+1,L) ) then

if(1) then
* write(*,*)'SOLUTION: setting plus neighbor less 1'
endif

neighbors( 2, P,L) = neighbors( 2, P,L) - 1

elseif( neighbors( 2, P,L) .lt.
. neighbors(1,P+1,L) ) then

if(1) then
* write(*,*)'SOLUTION: setting minus neighbor less 1'
endif

neighbors(1,P+1,L) = neighbors(1,P+1,L) - 1

else
* write(*,*)'ERROR: serious! layer,neibhors: ',L,
* . neighbors( 2, P,L),
* . neighbors(1,P+1,L)
endif

endif
enddo
 

check_sum = 0
do P = 1,peaks(L)
check_sum = check_sum + 1 + neighbors(1,P,L)+
. neighbors( 2,P,L)
enddo
if(check_sum .ne. hits(L) ) then
write(*,*) 'ERROR: layer,hits,check_sum: ',
. L,hits(L),check_sum

* dump = .true.

endif

enddo

* End of finding clusters...

* Find error weighted CoGs:

do L = 1,4
do P = 1,peaks(L)
 

sum_Q = 0.

do N = -neighbors(1,P,L),neighbors(2,P,L)

sum_Q = sum_Q + float(signal( peak_strip(P,L)+N,L))

enddo
 
 
 

if( 1+neighbors(1,P,L)+neighbors(2,P,L) .eq. 1) then

CoG(P,L) = float(strip(peak_strip(P,L),L))
e_CoG(P,L) = 1/sqrt(12.)

elseif(1+neighbors(1,P,L)+neighbors(2,P,L) .ge. 2) then

sum_N = 0.
sum_D = 0.

do N = -neighbors(1,P,L),neighbors(2,P,L)

sum_N = sum_N + float(strip( peak_strip(P,L)+N,L))*
. float(signal( peak_strip(P,L)+N,L))/
. float(noise( peak_strip(P,L)+N,L)**2)

sum_D = sum_D + float(signal( peak_strip(P,L)+N,L))/
. float(noise( peak_strip(P,L)+N,L)**2)

enddo

if( sum_D .ne. 0.) then
CoG(P,L) = sum_N / sum_D
else
* write(6,*) 'ERROR: sum_D .eq. 0.! n,d:',sum_n,sum_d
goto 999
endif
 

sum_E = 0.

do N = -neighbors(1,P,L),neighbors(2,P,L)

sum_E = sum_E + (
. ( float(strip( peak_strip(P,L)+N,L))-
. Cog(P,L)
. )/float(noise( peak_strip(P,L)+N,L))
. )**2

enddo

e_Cog(P,L) = sqrt(abs(sum_E))/abs(sum_D)
 
 

if(1) then

if( CoG(P,L) .lt.
. strip( peak_strip(P,L)-neighbors(1,P,L),L) .or.
. CoG(P,L) .gt.
. strip( peak_strip(P,L)+neighbors( 2,P,L),L) ) then

* write(*,*)'Event: ',event
* write(*,*) 'ERROR: bad COG, COG,strip:',CoG(P,L),
* . strip( peak_strip(P,L),L),
* . 1+neighbors(1,P,L)+neighbors(2,P,L),
* . e_CoG(P,L),sum_Q
goto 999
endif
endif
 
 
 
 
 

else
write(*,*) 'ERROR: illegal cluster hit sum: ',
. 1+neighbors(1,P,L)+neighbors(2,P,L)
goto 999
endif
enddo
enddo
 

* Now we are going to check for the copied hit problem::

copied_hit_candidates = 0

do L = 1,4
do P = 1,peaks(L)

copied(P,L) = .false.

enddo
enddo
 

do L1 = 1,3
do L2 = L1+1,4

do P1 = 1,peaks(L1)
do P2 = 1,peaks(L2)

if( neighbors( 2,P1,L1) .eq.
. neighbors( 2,P2,L2) .and.
. neighbors(1,P1,L1) .eq.
. neighbors(1,P2,L2) .and.
. strip(peak_strip(P1,L1),L1) .eq.
. strip(peak_strip(P2,L2),L2) ) then

copied(P1,L1) = .true.
copied(P2,L2) = .true.

copied_hit_candidates = copied_hit_candidates + 1
copied_hits = .true.

endif

enddo
enddo

enddo
enddo

* Find the number of 4 hit layer tracks...

n_trk_4(D) = 0

do P1 = 1,peaks(1)
do P2 = 1,peaks(2)
do P3 = 1,peaks(3)
do P4 = 1,peaks(4)

if( copied(P1,1) .or.
. copied(P2,2) .or.
. copied(P3,3) .or.
. copied(P4,4) ) then
goto 777
endif

n_trk_4(D) = n_trk_4(D) + 1

if( n_trk_4(D) .gt. 16 ) then
* write(*,*) 'ERROR: increase max_allowed_tracks: ',
* . n_trk_4(D),
* . peaks(1)*peaks(2)*peaks(3)*peaks(4),D
goto 999
endif
 

trk(1,n_trk_4(D),D) = CoG(P1,1) - offset(1,D)
trk(2,n_trk_4(D),D) = CoG(P2,2) - offset(2,D)
trk(3,n_trk_4(D),D) = CoG(P3,3) - offset(3,D)
trk(4,n_trk_4(D),D) = CoG(P4,4) - offset(4,D)

e_trk(1,n_trk_4(D),D) = e_CoG(P1,1)
e_trk(2,n_trk_4(D),D) = e_CoG(P2,2)
e_trk(3,n_trk_4(D),D) = e_CoG(P3,3)
e_trk(4,n_trk_4(D),D) = e_CoG(P4,4)

trk_peak(1,n_trk_4(D),D) = P1
trk_peak(2,n_trk_4(D),D) = P2
trk_peak(3,n_trk_4(D),D) = P3
trk_peak(4,n_trk_4(D),D) = P4

trk_layer(1,n_trk_4(D),D) = 1
trk_layer(2,n_trk_4(D),D) = 2
trk_layer(3,n_trk_4(D),D) = 3
trk_layer(4,n_trk_4(D),D) = 4
 

delta_12 = trk(1,n_trk_4(D),D) - trk(2,n_trk_4(D),D)
delta_23 = trk(2,n_trk_4(D),D) - trk(3,n_trk_4(D),D)
delta_34 = trk(3,n_trk_4(D),D) - trk(4,n_trk_4(D),D)
 

if(abs(delta_12-delta_23) .gt. max_two_layer_diff .or.
. abs(delta_23-delta_34) .gt. max_two_layer_diff ) then

n_trk_4(D) = n_trk_4(D) - 1

else

* Here is a measure of the "goodness" of the track (before the tilt correction).

delta_delta = abs(delta_12-delta_23) + abs(delta_23-delta_34)

do P1c = 1,peaks(1)
do P2c = 1,peaks(2)
do P3c = 1,peaks(3)
do P4c = 1,peaks(4)

if( copied(P1c,1) .or.
. copied(P2c,2) .or.
. copied(P3c,3) .or.
. copied(P4c,4) ) then
goto 888
endif

if( (P1c .eq. P1 .or.
. P2c .eq. P2 .or.
. P3c .eq. P3 .or.
. P4c .eq. P4 ) .and.
. .not.(P1c .eq. P1 .and.
. P2c .eq. P2 .and.
. P3c .eq. P3 .and.
. P4c .eq. P4 ) ) then
 

trk_c(1) = CoG(P1c,1) - offset(1,D)
trk_c(2) = CoG(P2c,2) - offset(2,D)
trk_c(3) = CoG(P3c,3) - offset(3,D)
trk_c(4) = CoG(P4c,4) - offset(4,D)
 

delta_12_check = trk_c(1) - trk_c(2)
delta_23_check = trk_c(2) - trk_c(3)
delta_34_check = trk_c(3) - trk_c(4)

if(abs(delta_12_check-delta_23_check) .le.
. max_two_layer_diff .and.
. abs(delta_23_check-delta_34_check) .le.
. max_two_layer_diff ) then
 

delta_delta_check = abs(delta_12_check-delta_23_check) +
. abs(delta_23_check-delta_34_check)
 

endif
endif
888 continue
enddo
enddo
enddo
enddo
 

endif
777 continue
enddo
enddo
enddo
enddo
 
 

if(n_trk_4(D) .gt. 1) then
* write(6,*) 'NOTE: more than 1 track,event,d:',
* . n_trk_4(D),event,d
endif

do T1 = 1,n_trk_4(D)-1
do T2 = T1+1,n_trk_4(D)

if( trk_peak(1,T1,D) .eq. trk_peak(1,T2,D) .or.
. trk_peak(2,T1,D) .eq. trk_peak(2,T2,D) .or.
. trk_peak(3,T1,D) .eq. trk_peak(3,T2,D) .or.
. trk_peak(4,T1,D) .eq. trk_peak(4,T2,D) ) then

* write(6,'(a,2i5,t60,8i3)')
* . 'WARNING: cluster conflict between tracks, evt,d',
* . event,d,trk_peak(1,T1,D), trk_peak(1,T2,D),
* . trk_peak(2,T1,D), trk_peak(2,T2,D),
* . trk_peak(3,T1,D), trk_peak(3,T2,D),
* . trk_peak(4,T1,D), trk_peak(4,T2,D)

* dump = .true.
endif
enddo
enddo
 

* Loop over all combinations of 3 layers...when?
 

n_trk_3(D) = 0

do L1 = 1,2
do L2 = L1+1,3
do L3 = L2+1,4

do P1 = 1,peaks(L1)
do P2 = 1,peaks(L2)
do P3 = 1,peaks(L3)

if( copied(P1,L1) .or.
. copied(P2,L2) .or.
. copied(P3,L3) ) then
goto 666
endif

do N = 1,n_trk_4(D)
do L = 1,n_layers
if(L1.eq.trk_layer(L,N,D).and.P1.eq.trk_peak(L,N,D) .or.
. L2.eq.trk_layer(L,N,D).and.P2.eq.trk_peak(L,N,D) .or.
. L3.eq.trk_layer(L,N,D).and.P3.eq.trk_peak(L,N,D) ) then
goto 666
endif
enddo
enddo
 

n_trk_3(D) = n_trk_3(D) + 1

if( n_trk_3(D) .gt. 16 ) then
* write(*,*) 'ERROR: 3 increase max_allowed_tracks: ',
* . n_trk_3(D),
* . peaks(L1)*peaks(L2)*peaks(L3),event,D
goto 999
endif
 

trk3(1,n_trk_3(D),D) = CoG(P1,L1) - offset(L1,D)
trk3(2,n_trk_3(D),D) = CoG(P2,L2) - offset(L2,D)
trk3(3,n_trk_3(D),D) = CoG(P3,L3) - offset(L3,D)

e_trk3(1,n_trk_3(D),D) = e_CoG(P1,L1)
e_trk3(2,n_trk_3(D),D) = e_CoG(P2,L2)
e_trk3(3,n_trk_3(D),D) = e_CoG(P3,L3)

trk_peak3(1,n_trk_3(D),D) = P1
trk_peak3(2,n_trk_3(D),D) = P2
trk_peak3(3,n_trk_3(D),D) = P3

trk_layer3(1,n_trk_3(D),D) = L1
trk_layer3(2,n_trk_3(D),D) = L2
trk_layer3(3,n_trk_3(D),D) = L3
 

delta_12 = trk3(1,n_trk_3(D),D) - trk3(2,n_trk_3(D),D)
delta_23 = trk3(2,n_trk_3(D),D) - trk3(3,n_trk_3(D),D)

if(abs(delta_12-delta_23) .gt. 2 * max_two_layer_diff ) then

n_trk_3(D) = n_trk_3(D) - 1

endif

666 continue
enddo
enddo
enddo

enddo
enddo
enddo
 

enddo ! of strip direction loop..
 

* write(*,*) 'event: ',event
* write(*,*) 'n_trk_4: ',n_trk_4
* write(*,*) 'n_trk_3: ',n_trk_3

if( n_trk_4( 2 ) .eq. 1 .and.
. n_trk_4( 1 ) .eq. 1 ) then
* write(*,*) 'Four: ',event
endif

if( n_trk_3( 2 ) .eq. 1 .and.
. n_trk_3( 1 ) .eq. 1 .and.
. n_trk_4( 2 ) .eq. 0 .and.
. n_trk_4( 1 ) .eq. 0 ) then
* write(*,*) 'Three:',event
endif

if((n_trk_3( 2 ) .eq. 1 .and.
. n_trk_3( 1 ) .eq. 0 .and.
. n_trk_4( 2 ) .eq. 0 .and.
. n_trk_4( 1 ) .eq. 1 ).or.
. (n_trk_3( 2 ) .eq. 0 .and.
. n_trk_3( 1 ) .eq. 1 .and.
. n_trk_4( 2 ) .eq. 1 .and.
. n_trk_4( 1 ) .eq. 0 ) ) then
* write(*,*) 'Mixed:',event
endif
 

if( .not.( ( n_trk_4( 2 ) .eq. 1 .and.
. n_trk_4( 1 ) .eq. 1
. )
. .or.( n_trk_3( 2 ) .eq. 1 .and.
. n_trk_3( 1 ) .eq. 1 .and.
. n_trk_4( 2 ) .eq. 0 .and.
. n_trk_4( 1 ) .eq. 0
. )
. .or.((n_trk_3( 2 ) .eq. 1 .and.
. n_trk_3( 1 ) .eq. 0 .and.
. n_trk_4( 2 ) .eq. 0 .and.
. n_trk_4( 1 ) .eq. 1 ).or.
. (n_trk_3( 2 ) .eq. 0 .and.
. n_trk_3( 1 ) .eq. 1 .and.
. n_trk_4( 2 ) .eq. 1 .and.
. n_trk_4( 1 ) .eq. 0 )
. )
. )
.
. ) goto 999

* OK, since we know n_trk_"x" is now zero or one:
 

do D = 1,2
do N = 1,n_trk_4(D)

do L = 1,4
tk(L) = trk(L,N,D)
etk(L) = e_trk(L,N,D)
enddo

if( D .eq. 1) then

call linear_fit(x_V, tk, etk, 4,
.
. Slope_V,
. Intercept_V,
. Error_slope_V,
. Error_intercept_V,
. Linear_correlation_constant_V,
. Least_squares_by_DOF_V)

else
call linear_fit(x_H, tk, etk, 4,
.
. Slope_H,
. Intercept_H,
. Error_slope_H,
. Error_intercept_H,
. Linear_correlation_constant_H,
. Least_squares_by_DOF_H)

endif
enddo
enddo
 
 

do D = 1,2
if( n_trk_4(D) .eq. 0) then
do N = 1,n_trk_3(D)

do L = 1,4
found_layer(L) = .false.
enddo
do L = 1,3
tk(L) = trk3(L,N,D)
etk(L) = e_trk3(L,N,D)
found_layer( trk_layer3(L,N,D) ) = .true.
enddo
do L = 1,4
if(.not. found_layer(L)) then
H = L
endif
enddo
 
 

if( D .eq. 1) then

do L = 1,4
if( L .lt. H) then
x3_V(L) = x_V(L)
elseif(L .gt. H) then
x3_V(L-1) = x_V(L)
endif
enddo
 

call linear_fit(x3_V, tk, etk, 3,
.
. Slope_V,
. Intercept_V,
. Error_slope_V,
. Error_intercept_V,
. Linear_correlation_constant_V,
. Least_squares_by_DOF_V)
 
 

else

do L = 1,4
if( L .lt. H) then
x3_H(L) = x_H(L)
elseif(L .gt. H) then
x3_H(L-1) = x_H(L)
endif
enddo
 

call linear_fit(x3_H, tk, etk, 3,
.
. Slope_H,
. Intercept_H,
. Error_slope_H,
. Error_intercept_H,
. Linear_correlation_constant_H,
. Least_squares_by_DOF_H)

endif
 

enddo
endif
enddo
 
 

do D = 1,2
do N = 1,n_trk_4(D)

if(D .eq. 1) then
do L = 1,4
trk(L,N,D) = trk(L,N,D)+(slope_h*x_V(L)+intercept_h)*
. tilt(L,D)/1023.
enddo
else
do L = 1,4
trk(L,N,D) = trk(L,N,D)+(slope_v*x_H(L)+intercept_v)*
. tilt(L,D)/1023.
enddo
endif

enddo

do N = 1,n_trk_3(D)

if(D .eq. 1) then
do L = 1,3

trk3(L,N,D) = trk3(L,N,D)+(slope_h*x3_V(L)+intercept_h)*
. tilt(trk_layer3(L,N,D),D)/1023.
enddo
else
do L = 1,3

trk3(L,N,D) = trk3(L,N,D)+(slope_v*x3_H(L)+intercept_v)*
. tilt(trk_layer3(L,N,D),D)/1023.
enddo
endif

enddo
enddo
 
 

* Now fit again, but with the tilt correction:

do D = 1,2
do N = 1,n_trk_4(D)
do L = 1,4
tk(L) = trk(L,N,D)
etk(L) = e_trk(L,N,D)
enddo
 

if( D .eq. 1) then

call linear_fit(x_V, tk, etk, 4,
.
. Slope_V,
. Intercept_V,
. Error_slope_V,
. Error_intercept_V,
. Linear_correlation_constant_V,
. Least_squares_by_DOF_V)
 

if(1) then
* write(*,*) 'Resolution: ',
* . sqrt( (x_P2*Error_slope_V)**2 +
* . (Error_intercept_V)**2 )
endif
 

else

call linear_fit(x_H, tk, etk, 4,
.
. Slope_H,
. Intercept_H,
. Error_slope_H,
. Error_intercept_H,
. Linear_correlation_constant_H,
. Least_squares_by_DOF_H)
 

if(1) then
* write(*,*) 'Resolution: ',
* . sqrt( (x_P2*Error_slope_H)**2 +
* . (Error_intercept_H)**2 )
endif

endif
 

enddo
enddo
 
 
 
 

do D = 1,2
if( n_trk_4(D) .eq. 0) then
do N = 1,n_trk_3(D)

do L = 1,4
found_layer(L) = .false.
enddo
do L = 1,3
tk(L) = trk3(L,N,D)
etk(L) = e_trk3(L,N,D)
found_layer( trk_layer3(L,N,D) ) = .true.
enddo
do L = 1,4
if(.not. found_layer(L)) then
H = L
endif
enddo
 

if( D .eq. 1) then

do L = 1,4
if( L .lt. H) then
x3_V(L) = x_V(L)
elseif(L .gt. H) then
x3_V(L-1) = x_V(L)
endif
enddo

call linear_fit(x3_V, tk, etk, 3,
.
. Slope_V,
. Intercept_V,
. Error_slope_V,
. Error_intercept_V,
. Linear_correlation_constant_V,
. Least_squares_by_DOF_V)

if(1) then

* write(*,*) x_v
* write(*,*) x3_v

* write(*,*) 'Resolution: ',
* . sqrt( (x_P2*Error_slope_V)**2 +
* . (Error_intercept_V)**2 )
endif
 

else

do L = 1,4
if( L .lt. H) then
x3_H(L) = x_H(L)
elseif(L .gt. H) then
x3_H(L-1) = x_H(L)
endif
enddo

call linear_fit(x3_H, tk, etk, 3,
.
. Slope_H,
. Intercept_H,
. Error_slope_H,
. Error_intercept_H,
. Linear_correlation_constant_H,
. Least_squares_by_DOF_H)

if(1) then

* write(*,*) x_h
* write(*,*) x3_h

*flag* write(*,*) 'Resolution: ',
* . sqrt( (x_P2*Error_slope_H)**2 +
* . (Error_intercept_H)**2 )

endif
endif
 

enddo
endif
enddo

* so, need to write out tracks to an array
 
 

do N = 1,2*4
final_t(N) = 9999.
e_final_t(N) = 9999.
enddo
 

do D = 1,2
do N = 1,n_trk_4(D)
do L = 1,4
final_t(4*(D-1)+L) = trk(L,N,D)
e_final_t(4*(D-1)+L) = e_trk(L,N,D)
enddo
enddo
enddo

do D = 1,2
if( n_trk_4(D) .eq. 0) then
do N = 1,n_trk_3(D)
do L = 1,4-1

final_t(4*(D-1)+ trk_layer3(L,N,D) ) = trk3(L,N,D)
e_final_t(4*(D-1)+ trk_layer3(L,N,D) ) = e_trk3(L,N,D)

enddo
enddo
endif
enddo
 

********************** Track Variables ********************

* . n_trk_4( 1),
* . n_trk_4(2),
* . n_trk_3( 1),
* . n_trk_3(2),
* . Slope_V,
* . Intercept_V,
* . Error_slope_V,
* . Error_intercept_V,
* . Linear_correlation_constant_V,
* . Least_squares_by_DOF_V,
* . Slope_H,
* . Intercept_H,
* . Error_slope_H,
* . Error_intercept_H,
* . Linear_correlation_constant_H,
* . Least_squares_by_DOF_H,
* . (final_t(N),N=1,8),
* . (e_final_t(N),N=1,8)
 

************************** Track Projection *****************
 

* Translates intercept into cm

Intercept_H = Intercept_H*0.005632
Intercept_V = Intercept_V*0.005632

Slope_H = Slope_Coefficient*Slope_H

* Extends tracks until it hits 3rd layer of chamber
Projection_H = Intercept_H + Slope_H*Travel_dist
Projection_V = Intercept_V + Slope_Coefficient*Slope_V*Travel_dist

* Pattern correction
open(79, FILE='pattern.dat')

do D = 1,255
read(79,*, end=300, err=200) pattern, temp_offset
if (pattern.eq.C_LCT_pattern) then
strip_offset = temp_offset
goto 300
endif

100 enddo
goto 300
200 write(*,*) 'ERROR!!!!'

300 close(79)

trans_mom = 1

if (c_lct_pattern.lt.192.and.c_lct_pattern.gt.173.or.
. c_lct_pattern.lt.142) then
strip_offset = 4*strip_offset
trans_mom = 0
endif
 
 

* Translates Projected value into the half strip it hits
if (trans_mom.eq.1) then
Half_strip = (Projection_H - C_LCT_offset)/Half_strip_width
. - strip_offset
else
Half_strip = ((Projection_H - C_LCT_offset)/
. Half_strip_width)/4. - 0.375
. - strip_offset*4.
endif

* Translates projected value into the anode wire it hits
Anode_wire = (Projection_V - A_LCT_offset)/Wire_spacing

* Subtracts the half strip hit on the CLCT board from the predicted half strip
HS_diff = Half_Strip - C_LCT_half_strip

* Same for anodes
Anode_diff = Anode_wire - A_LCT_wire
 

* Puts info into ntuple 20
ntuple(1) = event
ntuple(2) = Slope_H
ntuple(3) = Intercept_H
ntuple(4) = Slope_V
ntuple(5) = Intercept_V
ntuple(6) = Projection_H
ntuple(7) = Projection_V
ntuple(8) = Half_strip
ntuple(9) = Anode_wire
ntuple(10) = c_lct_half_strip
ntuple(11) = a_lct_wire
ntuple(12) = HS_diff
ntuple(13) = Anode_diff
ntuple(14) = trans_mom
ntuple(15) = c_lct_pattern

call hfn(20, ntuple)
 

999 continue

* OK = .true.

111 continue
return
end

SUBROUTINE LINEAR_FIT( X,
. Y,
. eY,
. Num_points,
. Slope,
. Intercept,
. Error_slope,
. Error_intercept,
. Linear_correlation_constant,
. Least_squares_by_DOF)
********************************************************************
* NOTE: each point is weighted by its error.
*
* NOTE: Y(i) = slope * X(i) + intercept
*
********************************************************************
IMPLICIT NONE

INTEGER
. Num_points

REAL
. X(*),
. Y(*),
. eY(*),
. slope,
. intercept,
. error_slope,
. error_intercept,
. linear_correlation_constant,
. Least_squares_by_DOF

REAL
. denominator,
. delta_x,
. delta_y,
. correlation,
. variance,
. sum_x,
. sum_y,
. sum_xx,
. sum_yy,
. sum_xy,
. sum,
. weight

INTEGER
. I

slope = 0.
intercept = 0.
error_slope = 0.
error_intercept = 0.
linear_correlation_constant = 0.
Least_squares_by_DOF = 0.

IF(Num_points .LT. 2) Then
* WRITE(*,*)'ERROR: LINEAR_FIT, not enough points!'
GOTO 999
Endif

do i = 1,Num_points
if(ey(i) .eq. 0.) then
* write(6,*) 'ERROR: weight is zero, so can''t do a linear fit!'
goto 999
endif
enddo
 

sum = 0.
sum_x = 0.
sum_y = 0.
sum_xx = 0.
sum_yy = 0.
sum_xy = 0.
 

DO i = 1,Num_points

weight = 1./ey(i)**2

sum = sum + weight
sum_x = sum_x + weight*x(i)
sum_y = sum_y + weight*y(i)
sum_xx = sum_xx + weight*x(i)**2
sum_yy = sum_yy + weight*y(i)**2
sum_xy = sum_xy + weight*x(i)*y(i)

ENDDO
 

delta_x = sum * sum_xx - sum_x * sum_x
delta_y = sum * sum_yy - sum_y * sum_y
correlation = sum * sum_xy - sum_x * sum_y
 

IF(delta_x .eq. 0.) Then
* WRITE(*,*)'ERROR: LINEAR_FIT, delta=0, infinite slope'
GOTO 999
Endif
 

intercept = (sum_xx * sum_y - sum_x * sum_xy)/delta_x

slope = correlation/delta_x
 

IF(Num_points .gt. 2) Then

variance = (sum_yy +
. sum*intercept**2 +
. sum_xx*slope**2 -
. 2*(intercept*sum_y +
. slope *sum_xy -
. slope*intercept*sum_x )) / (Num_points - 2)

IF(variance .lt. 0) Then
* WRITE(*,*)'ERROR:LINEAR_FIT: Changing sign of neg. variance:',
* . variance,
* . Num_points,slope,intercept,sum_x,sum_y,sum_xx,sum_yy,sum_xy
* WRITE(*,*)' Precision problems?'
variance = -variance
Endif

error_intercept = sqrt(variance*sum_xx/delta_x)
error_slope = sqrt(variance*sum/delta_x)
 

Least_squares_by_DOF = variance

Endif
 
 

denominator = delta_x*delta_y

IF(denominator .gt. 0.) Then

denominator = sqrt(denominator)
linear_correlation_constant = correlation/denominator

Else
IF(delta_y .ne. 0) Then
* WRITE(*,*)'ERROR: correlation denominator in LINEAR_FIT! ',
* . delta_x,delta_y,denominator,correlation
Endif
Endif
 

999 CONTINUE
RETURN
END
 
 

999 return
dump_tele = 1.
end
Quit
 
 

cd //res

nt/loop 10 dump_tele 51805 1

*type cd //pawc to access ntuple

RETURN