subroutine v_merge
C-----------------------------------------------------------------------
C-
C- Created 12-Jun-1996 Alexander Andryakov
C-
C- Purpose: Goes through all t.c. found by 1-view p.r. and chooses the pairs
C- to be merged in one 3-d t.c.
C-
C- Output :
C- numb_3d_cand - number of 3-d t.c. built
C----------------------------------------------------------------------
$$IMPLICIT
$$INCLUDE 'K$INC:BCS.INC'
$$INCLUDE 'YBOS$LIBRARY:ERRCOD.INC'
$$INCLUDE 'K$INC:ERLEVL.INC'
$$INCLUDE 'K$ITRK:CDTCUT.INC'
$$INCLUDE 'K$ITRK:DPRS.INC'
$$INCLUDE 'K$ITRK:DHRE.INC'
$$INCLUDE 'K$ITRK:DCMERGINT.INC'
c
INTEGER BLOCAT,BNUMB
LOGICAL intersect_in_layer_phi
c
INTEGER cycle,status,ind,
& inddat_dprs_1,itr1,inddat_dprs_2,itr2,view_n1,view_n2,
& checked_tc_list(N_TRK_MAX),itr_tried,
& merging_order(FPL_LNG_MAX),list_length
REAL first_pass_list(FPL_LNG_MAX,3),curvature1,curvature2
c
c
c Some preparations
c Take the value of LAYER_GAP and COMBINING_GAP (wire/phi gap)
c equal to corresponding cuts used in DTPTRN (one view pattern recognition)
c
LAYER_GAP = LDTXTR * 2
COMBINING_GAP = DTXPN3 + 0.01
c
c Locate DHRE bank
c
status = blocat(IW,'DHRE',1,ind,inddat_dhre)
if (status .ne. YESUCC) return ! DHRE bank has not been located
dhre_e_size = IW(inddat_dhre + DHRNCO)
inddat_dhre = inddat_dhre + IW(inddat_dhre + DHRHDS)
c
c FIRST loop over all t.c. Only t.c. from different stereo views and the
c same sign of curvatures are merged in accordance with ordered
c first_pass_list. If some t.c. becomes 3d one it leaves this loop
c and will be treated in the SECOND loop. (It's very important in
c multitrack regime) Start with the first t.c. in DPRS chain (take
c it as current)
c Build the ordered first_pass_list.
c
call make_1st_pass_list(first_pass_list,list_length,
& merging_order)
c
c Go through the list
c
if(list_length.lt.1) return
call vzero(checked_tc_list,N_TRK_MAX)
n_itr1_checked = 0
n_itr2_checked = 0
cycle=1
do while(cycle .le. list_length .and.
& n_itr1_checked .lt. first_itr2 - 1 .and.
& n_itr2_checked .lt. last_itr2-first_itr2+1)
c
c Get the t.c. attributes
c
itr1=int(first_pass_list(merging_order(cycle),2))
itr2=int(first_pass_list(merging_order(cycle),3))
inddat_dprs_1 = inddat_dprs_array(itr1)
inddat_dprs_2 = inddat_dprs_array(itr2)
c
c Check if neither of t.c. has already become 3d one
c
view_n1 = iw(inddat_dprs_1+DPRPOK)
view_n2 = iw(inddat_dprs_2+DPRPOK)
if(view_n1 .ne. 3 .and. view_n2 .ne. 3) then
c
c Try to merge them
c
if(inddat_dprs_1 .gt. 0 .and. inddat_dprs_2 .gt. 0) then
call merge_2_tc(itr1,inddat_dprs_1,itr2,inddat_dprs_2)
c
c Save for "curretnt" t.c. its partner t.c. number. It's done with the
c aim not to try to merge this pair in the SECOND loop
c
if(itr1 .le. N_TRK_MAX) checked_tc_list(itr1)=itr2
endif
endif
c
c Take the next pair of t.c.
c
cycle = cycle+1
enddo
c
c SECOND loop over all t.c. Only t.c. from different stereo views with the
c same sign and closest values of curvatures and NOT tried in the FIRST loop
c are about to be merged
c Locate the first t.c. in DPRS chain
c
status=BLOCAT(iw,'DPRS',-1,ind,inddat_dprs_1)
c
c Get the track candidate number
c
status=BNUMB(IW,ind,itr1)
if (status.ne.YESUCC) then
call ERLOGR('v_merge',ERSEVR,0,status,
& 'track number cannot be retrieved!?')
return
endif
c
c Go through all t.c.. (Loop over "current" t.c.)
c
do while(itr1 .lt. last_itr2)
if(inddat_dprs_array(itr1) .gt. 0) then
c
c DPRS bank of the "current" t.c. exists. Check all the rest t.c.
c for merging with it.
c
inddat_dprs_1=inddat_dprs_array(itr1)
c
c Get curvature and view number of the "current" t.c.
c
curvature1= rw(inddat_dprs_1+DPRCUR)
view_n1 = iw(inddat_dprs_1+DPRPOK)
c
c Get the number of t.c. already tried to be merged with the "current"
c one in the FIRST loop
c
if(itr1.le.N_TRK_MAX) then
itr_tried=checked_tc_list(itr1)
else
itr_tried=-1
endif
c
c Go throuhg all t.c. starting from the next to the "current" one and
c check them one by one if they have appropriate sign of curvature and
c view number (The last means that we will try to merge here only t.c.
c from different views 1,2,3. 3 means 3-d t.c.)
c Loop over "second" t.c..
c
itr2=itr1+1
do while(itr2 .le. last_itr2 .and.
& inddat_dprs_1 .gt. 0)
if(inddat_dprs_array(itr2) .gt. 0) then
c
c DPRS bank of "second" t.c. exists. Check it on merging
c
inddat_dprs_2=inddat_dprs_array(itr2)
c
c Get curvature and view number of "second" t.c.
c
curvature2= rw(inddat_dprs_2+DPRCUR)
view_n2 = iw(inddat_dprs_2+DPRPOK)
c
c Check if this pair of t.c. has not been already tried in the FIRST loop.
c This pair should not be attempted only if NEITHER of both t.c. has
c become 3d t.c. during the FIRST loop
c
if(itr2 .ne. itr_tried .or.
& view_n1 .eq. 3 .or. view_n2 .eq. 3) then
c
c Check whether these two t.c. are from different views
c and
c whether their curvatures are of the same sign
c and
c whether they occupy intersecting regions in layer-phi space
c and make an attempt to merge this pair of t.c..
c
c
if(curvature1 * curvature2 .gt. 0. .and.
& view_n1 .ne. view_n2 .and.
& intersect_in_layer_phi(inddat_dprs_1,inddat_dprs_2)
& ) then
call merge_2_tc(itr1,inddat_dprs_1,itr2
& ,inddat_dprs_2)
c
c Update curvature and view number of the "current" t.c.
c
inddat_dprs_1=inddat_dprs_array(itr1)
if(inddat_dprs_1 .gt. 0) then
curvature1= rw(inddat_dprs_1+DPRCUR)
view_n1 = iw(inddat_dprs_1+DPRPOK)
endif
endif
endif
endif
c
c Take next "second" t.c. and continue the loop over "second" t.c.
c
itr2 = itr2 + 1
enddo
c
c All the t.c. next to the "current" one have been checked/merged
c with it. Take next "current" t.c. and continue the loop over
c "current" t.c.
c
endif
itr1 = itr1 + 1
enddo
return
end
C===========================================================================
integer function adjust_layer_directions(
& nhits1,inddat_dprs_1,
& hitn_1_c_f,hitn_1_c_l,layer_1,layer_dir_1,
& nhits2,inddat_dprs_2,
& hitn_2_c_f,hitn_2_c_l,layer_2,layer_dir_2)
$$IMPLICIT
integer
& nhits1,inddat_dprs_1,
& hitn_1_c_f,hitn_1_c_l,layer_1,layer_dir_1,
& nhits2,inddat_dprs_2,
& hitn_2_c_f,hitn_2_c_l,layer_2,layer_dir_2
C-----------------------------------------------------------------------
C-
C- Created 11-Jun-1996 Alexander Andryakov
C-
C- Purpose: Finds the pieces of t.c. with the same layer direction of
C- movement.
C- It's supposed that at least one t.c. changes its direction
C-
C- Input :
C- nhits1 \ total number of hits of both t.c.
C- nhits2 /
C- inddat_dprs_1 \ pointer to the beginning of the hit list in
C- inddat_dprs_2 / DPRS banks of 1st and 2nd t.c. respectevly
C- hitn_1_c_f \ DHRE numbers of the 1st and the last hits in the
C- hitn_1_c_l / "current" cluster of the 1st t.c.
C- hitn_2_c_f \ same as hitn_1_c_f,hitn_1_c_l but for 2nd t.c.
C- hitn_2_c_l /
C- layer_1 \ layer number of clusters of both t.c.
C- layer_2 /
C- layer_dir_1 \ layer direction of both t.c.
C- layer_dir_2 /
C- Output :
C- hitn_1_c_f \
C- hitn_1_c_l | These values define layer cluster.
C- hitn_2_c_f | ONLY for one of 2 t.c. they will describe
C- hitn_2_c_l | new "current" cluster
C- layer_1 |
C- layer_2 /
C- Return :
C- <0 - adjustment failed
C- >0 - ---//--- succeeded
C----------------------------------------------------------------------
$$INCLUDE 'K$INC:BCS.INC'
$$INCLUDE 'YBOS$LIBRARY:ERRCOD.INC'
$$INCLUDE 'K$ITRK:DPRS.INC'
$$INCLUDE 'K$ITRK:DHRE.INC'
$$INCLUDE 'K$ITRK:DCMERGINT.INC'
c
INTEGER get_layer_cluster
INTEGER inddat_dprs,hitn_c_f,hitn_c_l,hitn_n_f,hitn_n_l,
& layer_dir,nhits,layer,
& layer_n
LOGICAL first
c
adjust_layer_directions = -1
if(iw(inddat_dprs_1 - DPRCEN + DPRLRG) .lt. 0 .and.
& iw(inddat_dprs_2 - DPRCEN + DPRLRG) .lt. 0) then
c
c Both t.c. change their direction
c For a moment we are only looking into the longest t.c. and are using only
c one layer direction matching and ignore another matching
c
else if(iw(inddat_dprs_1 - DPRCEN + DPRLRG) .lt. 0) then
if(nhits1 .ge. nhits2) then
first=.TRUE.
inddat_dprs = inddat_dprs_1
hitn_n_f = hitn_1_c_l+1
hitn_c_f = hitn_1_c_f
layer = layer_1
layer_dir = layer_dir_1
nhits = nhits1
else
first=.FALSE.
inddat_dprs = inddat_dprs_2
hitn_n_f = hitn_2_c_l+1
hitn_c_f = hitn_2_c_f
layer = layer_2
layer_dir = layer_dir_2
nhits = nhits2
endif
c
c Only 1st t.c. changes its direction.
c Find the beginning of its piece of other direction
c
first=.TRUE.
inddat_dprs = inddat_dprs_1
hitn_n_f = hitn_1_c_l+1
hitn_c_f = hitn_1_c_f
layer = layer_1
layer_dir = layer_dir_1
nhits = nhits1
else if(iw(inddat_dprs_2 - DPRCEN + DPRLRG) .lt. 0) then
c
c Only 2nd t.c. changes its direction
c Find the beginning of its piece of other direction
c
first=.FALSE.
inddat_dprs = inddat_dprs_2
hitn_n_f = hitn_2_c_l+1
hitn_c_f = hitn_2_c_f
layer = layer_2
layer_dir = layer_dir_2
nhits = nhits2
endif
c
c Loop over wires/clusters to find the piece of another direction
c
do while (hitn_n_f .le. nhits)
c
c Take next cluster and check new layer direction between this ckuster and
c previous one
c
hitn_n_l = get_layer_cluster(nhits,hitn_n_f,
& inddat_dprs,layer_n)
if((layer_n-layer)*layer_dir .lt. 0) then
adjust_layer_directions = 1
go to 111
endif
hitn_c_f = hitn_n_f
hitn_c_l = hitn_n_l
layer = layer_n
hitn_n_f = hitn_c_l + 1
enddo
111 if(first) then
hitn_1_c_f = hitn_c_f
hitn_1_c_l = hitn_c_l
layer_1 = layer
else
hitn_2_c_f = hitn_c_f
hitn_2_c_l = hitn_c_l
layer_2 = layer
endif
return
end
C===========================================================================
logical function check_azim_gap(
& hit1, inddat_dprs1, hit2, inddat_dprs2)
$$IMPLICIT
integer hit1, inddat_dprs1, hit2, inddat_dprs2
C----------------------------------------------------------------------
C-
C- Created 23-May-1996 Alexander Andryakov
C-
C- Input :
C- hit1 - hit number in DPRS bank of the 1st t.c.
C- hit2 - hit number in DPRS bank of the 2nd t.c.
C- inddat_dprs1 - pointer to the beginning of the data block of DPRS
C- bank of the 1st t.c.
C- inddat_dprs2 - the same as inddat_dprs1 but for the 2nd t.c.
C- Return :
C- result of azimuthal gap check (azimuthal gap is equivalent to the
C- gap measured in number of wires)
C----------------------------------------------------------------------
$$INCLUDE 'K$INC:BCS.INC'
$$INCLUDE 'YBOS$LIBRARY:ERRCOD.INC'
$$INCLUDE 'K$ITRK:DPRS.INC'
$$INCLUDE 'K$ITRK:DHRE.INC'
$$INCLUDE 'K$ITRK:CTIDRI.INC'
$$INCLUDE 'K$ITRK:CDTCUT.INC'
$$INCLUDE 'K$ITRK:CDTHIT.INC'
$$INCLUDE 'K$ITRK:CCONST.INC'
$$INCLUDE 'K$ITRK:DCMERGINT.INC'
c
INTEGER dhre_w_n1,dhre_w_n2,dhre_pos,layer1,layer2
REAL tan_stereo,azim_1,azim_2,azim_gap,phi1,phi2
c
dhre_w_n1 = iw(inddat_dprs1 + DPRLWI*(hit1-1)+DPRIWI)
dhre_pos = inddat_dhre + dhre_e_size*(dhre_w_n1-1)
layer1 = iw(dhre_pos + DHRSLR)
phi1 = TIDDP(layer1) * (iw(dhre_pos + DHRWNR)-1)
dhre_w_n2 = iw(inddat_dprs2 + DPRLWI*(hit2-1)+DPRIWI)
dhre_pos = inddat_dhre + dhre_e_size*(dhre_w_n2-1)
layer2 = iw(dhre_pos + DHRSLR)
phi2 = TIDDP(layer2) * (iw(dhre_pos + DHRWNR)-1)
if(mod(abs(layer1-layer2),2) .ne. 0) then
c
c The hits belong to layers of DIFFERENT views
c
tan_stereo=TIDST(layer1)/TIDSC(layer1)
azim_1=atan(tan_stereo/sqrt((TIDRA(layer1)/ZMAX)**2-
& tan_stereo**2))
tan_stereo=TIDST(layer2)/TIDSC(layer2)
azim_2=atan(tan_stereo/sqrt((TIDRA(layer2)/ZMAX)**2-
& tan_stereo**2))
azim_gap=abs(azim_1-azim_2)
else
c
c The hits belong to layers of the SAME view. Take the biggest azimuthal gap
c from one view pattern recognition.
c
azim_gap=max(TIDDP(layer1),TIDDP(layer2)) * DTXPN1
endif
if(abs(phi1-phi2).gt.azim_gap
& .and.
& pi2-abs(phi1-phi2).gt.azim_gap) then
check_azim_gap=.FALSE.
else
check_azim_gap=.TRUE.
endif
return
end
C===========================================================================
integer function combine_clusters(
& hitn_1_f,hitn_1_l,inddat_dprs_1,nhits1,
& hitn_2_f,hitn_2_l,inddat_dprs_2,nhits2,
& combine_buffer)
$$IMPLICIT
integer
& hitn_1_f,hitn_1_l,inddat_dprs_1,nhits1,
& hitn_2_f,hitn_2_l,inddat_dprs_2,nhits2,
& combine_buffer(1)
C----------------------------------------------------------------------
C-
C- Created 30-May-1996 Alexander Andryakov
C-
C- Input :
C- inddat_dprs_1 \ pointer to the beginning of the data block of DPRS
C- inddat_dprs_2 / banks of 1st and 2nd t.c. respectevly
C- hitn_1_f \ DHRE numbers of the 1st and the last hits in the
C- hitn_1_l / cluster of the 1st t.c.
C- hitn_2_f \ same as above but for 2nd t.c.
C- hitn_2_l /
C- Output :
C- combine_buffer - temporary buffer with mixture of wires from
C- both t.c.
C- Return :
C- <0 - combining failed
C- >0 - total number of wires saved into temporary buffer
C-
C- **** Very important note ***
C- It's supposed that both t.c. moves in the same azimuth direction
C- However the check of this is done
C----------------------------------------------------------------------
$$INCLUDE 'K$INC:BCS.INC'
$$INCLUDE 'YBOS$LIBRARY:ERRCOD.INC'
$$INCLUDE 'K$ITRK:CTIDRI.INC'
$$INCLUDE 'K$ITRK:DPRS.INC'
$$INCLUDE 'K$ITRK:DHRE.INC'
$$INCLUDE 'K$ITRK:DCMERGINT.INC'
c
INTEGER get_layer_cluster
LOGICAL check_azim_gap
c
INTEGER cycle1,cycle2, n_w, dhre_w_n, dhre_pos,layer,
& layer_1,hitn_1_n_f,hitn_1_n_l,
& layer_2,hitn_2_n_f,hitn_2_n_l
REAL phi_1_c, phi_1_l, phi_2_c, phi_2_l,
& phi_dir_1, phi_dir_2, phi_dir, azimthal_gap
c
REAL A,B,PROX
PROX(B,A)=B + 6.2831853 * ANINT(0.15915494*(A-B))
c
c
c Azimuth of 1st wire of 1st t.c. and delta phi of the layer
c
dhre_w_n = iw(inddat_dprs_1 + DPRLWI*(hitn_1_f-1)+DPRIWI)
dhre_pos = inddat_dhre + dhre_e_size*(dhre_w_n-1)
layer = iw(dhre_pos + DHRSLR)
phi_1_c = TIDDP(layer) * (iw(dhre_pos + DHRWNR)-1)
azimthal_gap=COMBINING_GAP*TIDDP(layer)
c
c Azimuth of last wire of 1st t.c.
c
dhre_w_n = iw(inddat_dprs_1 + DPRLWI*(hitn_1_l-1)+DPRIWI)
dhre_pos = inddat_dhre + dhre_e_size*(dhre_w_n-1)
phi_1_l = TIDDP(layer) * (iw(dhre_pos + DHRWNR)-1)
phi_dir_1 = sign(1.,phi_1_l-PROX(phi_1_c,phi_1_l))
c
c Azimuth of 1st wire of 2nd t.c.
c
dhre_w_n = iw(inddat_dprs_2 + DPRLWI*(hitn_2_f-1)+DPRIWI)
dhre_pos = inddat_dhre + dhre_e_size*(dhre_w_n-1)
phi_2_c = TIDDP(layer) * (iw(dhre_pos + DHRWNR)-1)
c
c Azimuth of last wire of 2nd t.c.
c
dhre_w_n = iw(inddat_dprs_2 + DPRLWI*(hitn_2_l-1)+DPRIWI)
dhre_pos = inddat_dhre + dhre_e_size*(dhre_w_n-1)
phi_2_l = TIDDP(layer) * (iw(dhre_pos + DHRWNR)-1)
phi_dir_2 = sign(1.,phi_2_l-PROX(phi_2_c,phi_2_l))
c
c Before performing combining procedure, one should check that this
c procedure will not destroy already existing t.c. It means that one should
c check azimuthal connection between the last wire of "current" cluster of
c one t.c. and the first wire of "next" cluster of other t.c. One can
c restrict oneself by checking only one such connection with the smallest
c layer gap.
c Take the "next" clusters of both t.c.
c
layer_1=0
hitn_1_n_f = hitn_1_l+1
if(hitn_1_n_f.le.nhits1)
& hitn_1_n_l = get_layer_cluster(
& nhits1,hitn_1_n_f,
& inddat_dprs_1,layer_1)
layer_2=0
hitn_2_n_f = hitn_2_l+1
if(hitn_2_n_f.le.nhits2)
& hitn_2_n_l = get_layer_cluster(
& nhits2,hitn_2_n_f,
& inddat_dprs_2,layer_2)
c
c 1st t.c. is 3d one and has "next" layer cluster.
c Check azimuth gap between hitn_2_l and hitn_1_n_f
c
if(layer_1 .gt. 0 .and. mod(layer_1-layer,2) .ne. 0) then
if(.not. check_azim_gap(
& hitn_1_n_f, inddat_dprs_1, hitn_2_l, inddat_dprs_2)) then
c
c Azimuthal gap is to large. Combining will damage 3d t.c.
c
combine_clusters = -1
go to 111
endif
else if(layer_2 .gt. 0 .and. mod(layer_2-layer,2) .ne. 0) then
c
c 2nd t.c. is 3d one and has "next" layer cluster.
c Check azimuth gap between hitn_1_l and hitn_2_n_f
c
if(.not. check_azim_gap(
& hitn_1_l, inddat_dprs_1, hitn_2_n_f, inddat_dprs_2)) then
c
c Azimuthal gap is to large. Combining will damage 3d t.c.
c
combine_clusters = -1
go to 111
endif
endif
c
c Now when we are sure that combining will not damage possible 3d t.c.
c let's start combining with checking azimuthal direction of movement and
c azimuthal overlap of combining clusters
c If at least one of the clusters has only 1 wire we do not need the
c following checks of movement directions coincidence, but only overlaping
c
if (hitn_1_l-hitn_1_f .gt. 0 .and.
& hitn_2_l-hitn_2_f .gt. 0) then
c
c Check the movement direction
c
if(phi_dir_1 .ne. phi_dir_2) then
c
c Different azimuth direction.
c Exit combining with unsuccessful return code
c
combine_clusters = -1
go to 111
endif
phi_dir = phi_dir_1
else if (hitn_1_l-hitn_1_f .eq. 0 .and.
& hitn_2_l-hitn_2_f .eq. 0) then
c
c Both t.c. has just single wires. Take azimuth direction of the movement
c from previous step.
c
if(phi_dir_cur_1 .ne. 0) then
phi_dir = phi_dir_cur_1
else if(phi_dir_cur_2 .ne. 0) then
phi_dir = phi_dir_cur_2
endif
else
c
c Only one t.c. has a single wire, but another has more than one. Define
c azimuth direction for this other t.c.
c
if (hitn_1_l-hitn_1_f .gt. 0) then
phi_dir = phi_dir_1
else
phi_dir = phi_dir_2
endif
endif
c
c Azimuth direction has been defined.
c Check now if t.c.'s overlap in azimuth
c
if(phi_dir*(phi_1_c-PROX(phi_2_l,phi_1_c)) .gt. azimthal_gap
& .or.
& phi_dir*(phi_2_c-PROX(phi_1_l,phi_2_c)) .gt. azimthal_gap)
& then
c
c No overlaping in azimuth.
c Exit combining with unsuccessful return code
c
combine_clusters = -1
go to 111
endif
c
c The same direction in azimuth & overlaping
c Proceed with combining, which at that point consists in putting ALL the
c wires of layer clusters of both t.c. into combining buffer in RIGHT ODER
c
cycle1 = hitn_1_f
cycle2 = hitn_2_f
n_w = 0
do while (cycle1.le.hitn_1_l.and.cycle2.le.hitn_2_l.and.
& n_w.lt.10000)
c
c We always start/continue at this point with current wires of both t.c.
c
dhre_w_n = iw(inddat_dprs_1 + DPRLWI*(cycle1-1)+DPRIWI)
dhre_pos = inddat_dhre + dhre_e_size*(dhre_w_n-1)
phi_1_c = TIDDP(layer) * (iw(dhre_pos + DHRWNR)-1)
dhre_w_n = iw(inddat_dprs_2 + DPRLWI*(cycle2-1)+DPRIWI)
dhre_pos = inddat_dhre + dhre_e_size*(dhre_w_n-1)
phi_2_c = TIDDP(layer) * (iw(dhre_pos + DHRWNR)-1)
c
c Now, which wire should we save into buffer and which t.c. should we take
c next wire from to proceed.
c
if(phi_dir*(phi_1_c-PROX(phi_2_c,phi_1_c)) .lt. 0) then
c
c phi_1_c < phi_2_c in a sence that it's before along the trajectory
c Save current wire and take next one from 1st t.c.
c
n_w = n_w + 1
combine_buffer(n_w)=
& iw(inddat_dprs_1 + DPRLWI*(cycle1-1)+DPRIWI)
cycle1 = cycle1 + 1
else
c
c phi_2_c < phi_1_c
c Save current wire and take next one from 2nd t.c.
c
n_w = n_w + 1
combine_buffer(n_w)=
& iw(inddat_dprs_2 + DPRLWI*(cycle2-1)+DPRIWI)
cycle2 = cycle2 + 1
endif
enddo
c
c The main loop is finished. There only two possibilites:
c 1. the buffer is full - nothing to do
c 2. one of t.c. clusters is over -
c saved the rest of another, starting from its current wire
c but only in case if combining succeeded (n_w>0)
c
if(n_w.eq.0) then
combine_clusters = -1
go to 111
endif
if(n_w.lt.10000) then
if(cycle1.le.hitn_1_l) then
c
c 2nd t.c. is over, save the rest of 1nd one
c
c********************************************** debugging only
if(cycle2.le.hitn_2_l) then
write (6,'('' after combinning ..... cycle1,cycle2,n_w:'',
& 3i10)') cycle1,cycle2,n_w
stop
endif
c********************************************** debugging only
do while (cycle1.le.hitn_1_l)
n_w = n_w + 1
combine_buffer(n_w)=
& iw(inddat_dprs_1 + DPRLWI*(cycle1-1)+DPRIWI)
cycle1 = cycle1 + 1
enddo
else if(cycle2.le.hitn_2_l) then
c
c 1st t.c. is over, save the rest of 2nd one
c
c********************************************** debugging only
if(cycle1.le.hitn_1_l) then
write (6,'('' after combinning ..... cycle1,cycle2,n_w:'',
& 3i10)') cycle1,cycle2,n_w
stop
endif
c********************************************** debugging only
do while (cycle2.le.hitn_2_l)
n_w = n_w + 1
combine_buffer(n_w)=
& iw(inddat_dprs_2 + DPRLWI*(cycle2-1)+DPRIWI)
cycle2 = cycle2 + 1
enddo
endif
endif
combine_clusters = n_w
111 return
end
C===========================================================================
integer function get_layer_cluster (
& nhits, cur_hit_n, inddat_dprs,
& cluster_layer)
$$IMPLICIT
integer nhits, cur_hit_n, inddat_dprs,
& cluster_layer
C----------------------------------------------------------------------
C-
C- Created 22-May-1996 Alexander Andryakov
C-
C- Input :
C- nhits - total number of hits of the t.c.
C- cur_hit_n - number of current hit of the t.c.
C- inddat_dprs - pointer to the beginning of the data block of DPRS
C- bank of the t.c.
C- Output:
C- cluster_layer - layer number of the cluster
C- Return :
C- number in DPRS of the last wire of found cluster
C----------------------------------------------------------------------
$$INCLUDE 'K$INC:BCS.INC'
$$INCLUDE 'YBOS$LIBRARY:ERRCOD.INC'
$$INCLUDE 'K$ITRK:DPRS.INC'
$$INCLUDE 'K$ITRK:DHRE.INC'
$$INCLUDE 'K$ITRK:DCMERGINT.INC'
c
INTEGER cycle,layer_c,dhre_w_n,dhre_pos,cluster_lenght
c
dhre_w_n = iw(inddat_dprs + DPRLWI*(cur_hit_n-1)+DPRIWI)
dhre_pos = inddat_dhre + dhre_e_size*(dhre_w_n-1)
cluster_layer = iw(dhre_pos + DHRSLR)
cluster_lenght=1
cycle = cur_hit_n
do while (cycle .lt. nhits)
cycle = cycle +1
dhre_w_n = iw(inddat_dprs + DPRLWI*(cycle-1)+DPRIWI)
dhre_pos = inddat_dhre + dhre_e_size*(dhre_w_n-1)
layer_c = iw(dhre_pos + DHRSLR)
if(layer_c .ne. cluster_layer) go to 111
cluster_lenght=cluster_lenght+1
enddo
111 get_layer_cluster=cur_hit_n+cluster_lenght-1
return
end
C===========================================================================
integer function get_layer_supercluster (
& nhits, cur_hit_n, inddat_dprs,bound_layer,layer_dir,
& cluster_layer,same_layer_f)
$$IMPLICIT
integer nhits, cur_hit_n, inddat_dprs,bound_layer,layer_dir,
& cluster_layer
logical same_layer_f
C----------------------------------------------------------------------
C-
C- Created 22-May-1996 Alexander Andryakov
C-
C- Input :
C- nhits - total number of hits of the t.c.
C- cur_hit_n - number of current hit of the t.c.
C- inddat_dprs - pointer to the beginning of the data block of DPRS
C- bank of the t.c.
C- bound_layer - layer number of the "next" cluster of other t.c.
C- layer_dir - layer direction of movement the t.c. piece
C- Output:
C- cluster_layer - layer number of the last cluster
C- same_layer_f - logical flag if the supercluster ends in the
C- bound_layer
C- Return :
C- number in DPRS of the last wire of found supercluster
C- <0 if no supercluster extension has been found
C----------------------------------------------------------------------
c
INTEGER get_layer_cluster
INTEGER next,int_last,int_layer
c
same_layer_f=.FALSE.
get_layer_supercluster=-1
next=cur_hit_n+1
do while (next .le. nhits)
int_last = get_layer_cluster(
& nhits,next,
& inddat_dprs,int_layer)
if(bound_layer. eq. int_layer) then
same_layer_f=.TRUE.
get_layer_supercluster=next-1
go to 111
else if(layer_dir*(bound_layer - int_layer) .gt. 0) then
cluster_layer = int_layer
get_layer_supercluster=int_last
next=int_last+1
else
go to 111
endif
enddo
111 return
end
C===========================================================================
integer function insert_inbetween(
& layer_1_c,layer_1_n,layer_2_c,
& inddat_dprs_1,nhits1,inddat_dprs_2,nhits2,
& hitn_1_c_f,hitn_1_c_l,hitn_1_n_f,hitn_1_n_l,
& hitn_2_c_f,hitn_2_c_l,
& num_w_saved,same_layer_f)
$$IMPLICIT
integer
& layer_1_c,layer_1_n,layer_2_c,
& inddat_dprs_1,nhits1,inddat_dprs_2,nhits2,
& hitn_1_c_f,hitn_1_c_l,hitn_1_n_f,hitn_1_n_l,
& hitn_2_c_f,hitn_2_c_l,num_w_saved
logical same_layer_f
C-----------------------------------------------------------------------
C-
C- Created 29-May-1996 Alexander Andryakov
C-
C- Purpose: Inserts the layer cluster (or supercluster) of 2nd
C- t.c. inbetween to clusters of 1st t.c.
C- *** Very important note : the numbering of t.c. is internal only
C-
C- Input :
C- layer_1_c - layer number of "current" cluster of 1st t.c.
C- layer_1_n - layer number of "next" cluster of 1st t.c.
C- layer_2_c - layer number of "current" cluster of 2nd t.c.
C- inddat_dprs_1 \ pointer to the beginning of the hit list of DPRS
C- inddat_dprs_2 / banks of 1st and 2nd t.c. respectevly
C- nhits1 - total number of hits of the 1st t.c.
C- nhits2 - total number of hits of the 2nd t.c.
C- hitn_1_c_f \ DHRE numbers of the 1st and the last hits in the
C- hitn_1_c_l / "current" cluster of the 1st t.c.
C- hitn_1_n_f \ the same as above but for the "next" cluster
C- hitn_1_n_l /
C- hitn_2_c_f \ same as hitn_1_c_f,hitn_1_c_l but for 2nd t.c.
C- hitn_2_c_l /
C- num_w_saved - total number of wires already saved into temporary
C- buffer of new merged t.c.
C- Output :
C- layer_2_c - layer number of last hit of 2nd t.c. supercluster
C- hitn_2_c_l - last hit number of 2nd t.c. supercluster
C- num_w_saved - total number of wires saved into temporary buffer
C- after the insertion
C- same_layer_f - logical flag (TRUE if the last clusters are in the
C- same layer
C- Return :
C- -1 insertion failed
C- 0 insertion has been done as required
C- 1 special case of successful insertion, when the last cluster of
C- 2nd t.c. supercluster and the "next" cluster of 1st t.c. are
C- in the SAME LAYER
C----------------------------------------------------------------------
$$INCLUDE 'K$ITRK:DCMERGINT.INC'
c
INTEGER get_layer_cluster,get_layer_supercluster,
& combine_clusters,save_clust2buf
LOGICAL check_azim_gap
c
INTEGER int_first,interim_last,int_layer,
& hitn_for_comb_l,l_f_comb,
& n_comb_wires,combine_buffer(BUFF_LENGTH),cycle
c
same_layer_f=.FALSE.
interim_last=-1
c
c First check what are the layer gaps between the current clusters of
c different t.c. (1st connection)
c
if(abs(layer_1_c - layer_2_c) .le. LAYER_GAP) then
c
c Layer gap of the first connection (1st cluster of 1st t.c. and
c cluster of 2nd t.c.) is ok.
c Check now azimuthal distance between last wire of 1st t.c. cluster
c and first wire of 2nd t.c. cluster Azimuthal gap is not constant
c and depends on the layer numbers.
c
if(check_azim_gap(hitn_1_c_l,inddat_dprs_1,
& hitn_2_c_f,inddat_dprs_2)) then
c
c Azimutal gap is ok. The 1st connection is established.
c Now try establish 2nd connection, one between the last wire of the
c cluster of the 2nd t.c. and the first wire of the 2nd cluster of
c the 1st t.c.
c First extend the current cluster of 2nd t.c. to supercluster
c
if(hitn_2_c_l+1.le.nhits2) then
int_first = hitn_2_c_l+1
int_layer = layer_2_c
interim_last = get_layer_supercluster(
& nhits2,hitn_2_c_l,
& inddat_dprs_2,layer_1_n,
& layer_1_n-layer_1_c,
& int_layer,same_layer_f)
endif
c
c Check if the last cluster of 2nd t.c. supercluster is in the same
c layer with the "next" cluster of 1st t.c.
c
if(same_layer_f) then
c
c Combine the hits of the same layer clusters. If they are combined
c successfully the insertion of the of the whole supercluster
c should be done.
c Remember that if same_layer_f=.TRUE. the get_layer_supercluster function
c returns results for supercluster WITHOUT this last cluster.
c So first get this cluster.
c
hitn_for_comb_l = get_layer_cluster(
& nhits2,interim_last+1,
& inddat_dprs_2,l_f_comb)
c
c Then perform combining
c
n_comb_wires=combine_clusters(
& hitn_1_n_f,hitn_1_n_l,inddat_dprs_1,nhits1,
& interim_last+1,hitn_for_comb_l,inddat_dprs_2,nhits2,
& combine_buffer)
if(n_comb_wires.gt.0) then
c
c The clusters of the same layer have been succsessfully combined.
c Save current clusters of 1st t.c., the whole supercluster of 2nd
c t.c. and and the second cluster of 1st t.c.
c 1.Save current clusters of 1st into temporary buffer
c
num_w_saved=save_clust2buf(
& inddat_dprs_1,hitn_1_c_f,hitn_1_c_l,
& num_w_saved)
c
c 2.Save the whole supercluster of 2nd up to the last "same" layer
c cluster
c
num_w_saved=save_clust2buf(
& inddat_dprs_2,hitn_2_c_f,interim_last,
& num_w_saved)
c
c 3.Save combined cluster into temporary buffer
c
cycle=1
do while(cycle.le.n_comb_wires .and.
& num_w_saved .lt. BUFF_LENGTH)
num_w_saved = num_w_saved +1
merge_buff(num_w_saved) = combine_buffer(cycle)
cycle = cycle + 1
enddo
c
c 4.Set successful return code for and exit
c
insert_inbetween=0
else
c Combining failed. Second connection is NOT established.
c Set unsuccessful return code and exit
c (
c It's probably to firm logic, but it's used here as the first approach.
c It can be loosened in the following way: if the above is not satisfied,
c one can instead try to connect the previous cluster from supercluster
c with the next cluster of partner t.c.. If the last is possible then the
c whole supercluster will be inserted inbetween two clusters and the last
c cluster of supercluster (that from "the same layer") will be dropped.
c )
c
insert_inbetween=-1
endif
c
c Redefine the upper boundary of supercluster of 2nd t.c.
c
if(interim_last.gt.0) then
hitn_2_c_l = hitn_for_comb_l
layer_2_c = l_f_comb
endif
* go to 111
else
c
c The supercluster is in between the two clusters.
c Redefine the upper boundary of the current cluster of 2nd t.c. to that
c of supercluster
c
if(interim_last.gt.0) then
hitn_2_c_l = interim_last
layer_2_c = int_layer
endif
c
c Check the second connection. Start with check of layer gap
c
if(abs(layer_2_c - layer_1_n) .le. LAYER_GAP) then
c
c Layer gap of the 2nd connection is ok.
c Check now azimuthal gap
c
if(check_azim_gap(hitn_2_c_l,inddat_dprs_2,
& hitn_1_n_f,inddat_dprs_1)) then
c
c Azimutal gap is ok. The 2nd connection is established.
c 1.Save current clusters of 1st & 2nd t.c. into temporary buffer
c
num_w_saved=save_clust2buf(
& inddat_dprs_1,hitn_1_c_f,hitn_1_c_l,
& num_w_saved)
num_w_saved=save_clust2buf(
& inddat_dprs_2,hitn_2_c_f,hitn_2_c_l,
& num_w_saved)
c
c 2.Set successful return code and exit
c
insert_inbetween=0
* go to 111
else
c
c Azimutal gap is too large. The 2nd connection is NOT established.
c Set unsuccessful return code and exit
c
insert_inbetween=-1
* go to 111
endif
else
c
c Layer gap is to large. The 2nd connection is NOT established.
c Set unsuccessful return code and exit
c
insert_inbetween=-1
* go to 111
endif
endif
else
c
c Azimutal gap is to large. The 1st connection is NOT established.
c Set unsuccessful return code and exit
c
insert_inbetween=-1
* go to 111
endif
else
c
c Layer gap is to large. The 1st connection is NOT established.
c Set unsuccessful return code and exit
c
insert_inbetween=-1
* go to 111
endif
111 return
end
C===========================================================================
logical function intersect_in_layer_phi(
& inddat_dprs1, inddat_dprs2)
$$IMPLICIT
integer inddat_dprs1, inddat_dprs2
C----------------------------------------------------------------------
C-
C- Created 24-Jun-1996 Alexander Andryakov
C-
C- Input :
C- inddat_dprs1 - pointer to the beginning of the data block of DPRS
C- bank of the 1st t.c.
C- inddat_dprs2 - the same as inddat_dprs1 but for the 2nd t.c.
C- Return :
C- result of check whether the regions of these 2 t.c. intersect
C----------------------------------------------------------------------
$$INCLUDE 'K$INC:BCS.INC'
$$INCLUDE 'YBOS$LIBRARY:ERRCOD.INC'
$$INCLUDE 'K$ITRK:CTIDRI.INC'
$$INCLUDE 'K$ITRK:CCONST.INC'
$$INCLUDE 'K$ITRK:DPRS.INC'
$$INCLUDE 'K$ITRK:DCMERGINT.INC'
c
INTEGER layer1_min,layer1_max,layer2_min,layer2_max,ltemp,
& layer_base, layer_baseu
REAL phi1_first,phi1_last,phi2_first,phi2_last,
& tan_stereo,azim_1,azim_2,azim_gap
c
REAL A,B,PROX
PROX(B,A)=B + PI2 * ANINT(0.15915494*(A-B))
c
intersect_in_layer_phi = .FALSE.
c
c Get layer boundaries of 1st t.c. region
c
ltemp = mod(abs(iw(inddat_dprs1+DPRLRG)),10000)
layer1_max = ltemp/100
layer1_min = mod(ltemp,100)
c
c Get layer boundaries of 2nd t.c. region
c
ltemp=mod(abs(iw(inddat_dprs2+DPRLRG)),10000)
layer2_max=ltemp/100
layer2_min=mod(ltemp,100)
c
c Check region "intersection" in layer dimension
c
if(layer1_min-layer2_max .le. LAYER_GAP .and.
& layer2_min-layer1_max .le. LAYER_GAP) then
c
c T.c. regions "intersect" in layer dimension
c Check if they "intersect" in azimuth dimension
c Get azimuthal boundaries of t.c.
c
phi1_first = rw(inddat_dprs1+DPRPHF)
phi1_last = rw(inddat_dprs1+DPRPHL)
phi2_first = rw(inddat_dprs2+DPRPHF)
phi2_last = rw(inddat_dprs2+DPRPHL)
if(phi1_last .lt. -100. .or. phi2_last .lt. -100.) then
c
c One of 2 t.c. occupy 2pi in azimuthal space
c
intersect_in_layer_phi = .TRUE.
else
c
c Calculate max allowed azimuthal gap
c Get minimal value of upper layer boundaries of t.c.
c
layer_base = min (layer1_max,layer2_max)
c
c Calculate max allowed azimuthal gap between layers:
c layer_base and layer_base+LAYER_GAP (if layer_base equal to NLAYER
c take layer_base-LAYER_GAP)
c
if(layer_base .lt. NLAYER) then
layer_baseu = layer_base + 1
else
layer_baseu = layer_base - 1
endif
tan_stereo=TIDST(layer_base)/TIDSC(layer_base)
azim_1=atan(tan_stereo/sqrt((TIDRA(layer_base)/ZMAX)**2-
& tan_stereo**2))
tan_stereo=TIDST(layer_baseu)/TIDSC(layer_baseu)
azim_2=atan(tan_stereo/sqrt((TIDRA(layer_baseu)/ZMAX)**2-
& tan_stereo**2))
azim_gap=abs(azim_1-azim_2)
c
c Now check azimuthal "intersection"
c
intersect_in_layer_phi =
& phi1_first - phi2_last .le. azim_gap .or.
& phi1_first - phi2_last - PI2 .le. azim_gap
& .or.
& phi2_first - phi1_last .le. azim_gap .or.
& phi2_first - phi1_last - PI2 .le. azim_gap
endif
endif
return
end
C===========================================================================
subroutine make_1st_pass_list(first_pass_list,list_length
& ,merging_order)
$$IMPLICIT
$$INCLUDE 'K$ITRK:DCMERGINT.INC'
real first_pass_list(FPL_LNG_MAX,3)
integer list_length,merging_order(FPL_LNG_MAX)
C---------------------------------------------------------------------------
C-
C- Created 7-Oct-1996 Alexander Andryakov
C-
C - Purpose: Goes through all t.c. found by 1-view p.r. and
C bulds the merging order list of t.c. pairs for the 1st
C merging loop
C-
C- Output :
C- first_pass_list - t.c. pairs list in the oreder of merging
C- (i,1) - "difference" order parameter
C- (i,2) - itr of 1st view t.c.
C- (i,3) - itr of 2nd view t.c.
C- list_length - number of t.c. combinations in the list
C- merging_order - array of pointers to first_pass_list elements
C- after ordering with SORTZV routine of CERNLIB
C---------------------------------------------------------------------------
$$INCLUDE 'K$INC:BCS.INC'
$$INCLUDE 'YBOS$LIBRARY:ERRCOD.INC'
$$INCLUDE 'K$INC:ERLEVL.INC'
$$INCLUDE 'K$ITRK:CDTCUT.INC'
$$INCLUDE 'K$ITRK:DPRS.INC'
$$INCLUDE 'K$ITRK:DHRE.INC'
c
LOGICAL intersect_in_layer_phi
c
INTEGER inddat_dprs_1,itr1,inddat_dprs_2,itr2,view_n1,view_n2
REAL difference,curvature1,curvature2,qual1,qual2
c
c LOOP over all t.c. Only t.c. from different stereo views and the
c same sign of curvatures are discussed.
c
call vzero(first_pass_list,3*FPL_LNG_MAX)
call vzero(merging_order,FPL_LNG_MAX)
list_length=0
itr1=1
do while(itr1 .lt. first_itr2 .and.
& inddat_dprs_array(itr1) .gt. 0 .and.
& list_length .lt. FPL_LNG_MAX)
inddat_dprs_1=inddat_dprs_array(itr1)
c
c Get curvature and view number of current t.c.
c
curvature1= rw(inddat_dprs_1+DPRCUR)
************************************************* debug only ***** beg
view_n1 = iw(inddat_dprs_1+DPRPOK)
if(view_n1 .ne. 1) then
write (6,'('' view_n1 = '',i10)') view_n1
stop
endif
************************************************* debug only ***** end
qual1 = rw(inddat_dprs_1+DPRQUA)
c
c Locate now all other t.c. of another view, check them one by one
c if they have appropriate sign of curvature and view number and see
c difference in curvatures
c
itr2 = first_itr2
do while(itr2 .le. last_itr2 .and.
& inddat_dprs_array(itr2) .gt. 0 .and.
& list_length .lt. FPL_LNG_MAX)
inddat_dprs_2=inddat_dprs_array(itr2)
c
c Get curvature and view number of next t.c.
c
curvature2= rw(inddat_dprs_2+DPRCUR)
************************************************* debug only ***** beg
view_n2 = iw(inddat_dprs_2+DPRPOK)
if(view_n2 .ne. 2) then
write (6,'('' view_n2 = '',i10)') view_n2
stop
endif
************************************************* debug only ***** end
c
c Check whether these two t.c. are from different views
c and
c whether their curvatures are of the same sign
c and
c whether they occupy intersecting regions in layer-phi space
c
if(curvature1 * curvature2 .gt. 0. .and.
& intersect_in_layer_phi(inddat_dprs_1,inddat_dprs_2)) then
c
c Build what serves as measure of how these t.c. are unlike to each
c other
c
qual2 = rw(inddat_dprs_2+DPRQUA)
c difference = abs(curvature1-curvature2) +
c & abs(qual1 - qual2)
difference = 2.*abs(curvature1-curvature2)/
& abs(curvature1+curvature2) +
& 2.*abs(qual1 - qual2)/abs(qual1 + qual2)
c
c Put these pair of t.c. with their "difference" value into list
c
list_length = list_length + 1
first_pass_list(list_length,1) = difference
first_pass_list(list_length,2) = float(itr1)
first_pass_list(list_length,3) = float(itr2)
endif
c
c Take next t.c. in the second stereo view
c
itr2 = itr2 + 1
enddo
c
c Take next t.c. in the first stereo view
c
itr1 = itr1 + 1
enddo
c
c Sort list in ascending order of "difference" value
c
call SORTZV(first_pass_list(1,1),merging_order,list_length,1,0,0)
return
end
C===========================================================================
subroutine merge_1_into_2(
& itr1,nhits1,inddat_dprs_1,
& hitn_1_c_f,hitn_1_c_l,layer_1_c,
& itr2,nhits2,inddat_dprs_2,
& hitn_2_c_f,hitn_2_c_l,layer_2_c,
& hitn_2_n_f,hitn_2_n_l,layer_2_n)
$$IMPLICIT
integer
& itr1,nhits1,inddat_dprs_1,
& hitn_1_c_f,hitn_1_c_l,layer_1_c,
& itr2,nhits2,inddat_dprs_2,
& hitn_2_c_f,hitn_2_c_l,layer_2_c,
& hitn_2_n_f,hitn_2_n_l,layer_2_n
C-----------------------------------------------------------------------
C-
C- Created 31-May-1996 Alexander Andryakov
C-
C- Purpose: Merge 1st t.c. into 2nd one.
C- 1st t.c. is always a single layer t.c.
C- Input :
C- itr1 \ DPRS bank nubmers of both t.c.
C- itr2 /
C- layer_1_c - layer number of "current" cluster of 1st t.c.
C- layer_2_c - layer number of "current" cluster of 2nd t.c.
C- layer_2_n - layer number of "next" cluster of 2nd t.c.
C- inddat_dprs_1 \ pointer to the beginning of the hit list of DPRS
C- inddat_dprs_2 / banks of 1st and 2nd t.c. respectevly
C- nhits1 \ total number of hits of
C- nhits2 / both t.c.
C- hitn_1_c_f \ DHRE numbers of the 1st and the last hits in the
C- hitn_1_c_l / "current" cluster of the 1st t.c.
C- hitn_2_c_f \ same as hitn_1_c_f,hitn_1_c_l but for 2nd t.c.
C- hitn_2_c_l /
C- hitn_2_n_f \ the same as above but for the "next" cluster
C- hitn_2_n_l /
C----------------------------------------------------------------------
$$INCLUDE 'K$INC:BCS.INC'
$$INCLUDE 'YBOS$LIBRARY:ERRCOD.INC'
$$INCLUDE 'K$ITRK:DPRS.INC'
$$INCLUDE 'K$ITRK:DHRE.INC'
$$INCLUDE 'K$ITRK:DCMERGINT.INC'
C
INTEGER get_layer_cluster,combine_clusters,
& insert_inbetween,modify_DPRS_DHRE
c
INTEGER
& num_w_saved,insert_status,layer_dir_2,
& n_comb_wires,combine_buffer(BUFF_LENGTH),cycle
LOGICAL same_layer_f
c
layer_dir_2 = layer_2_n - layer_2_c
num_w_saved=0
c
c Check if 1st t.c. should be really merged into the second one.
c
if(layer_dir_2*(layer_2_c-layer_1_c) .lt. 0) then
c
c 1st t.c. is AFTER "current" cluster of 2nd t.c.
c
c Loop over wires/clusters of 2nd t.c. (hit lists merging) to find two
c cluster between which 1st t.c. should be inserted
c
do while (hitn_2_c_f .le. nhits2)
c
c Check if the 1st t.c. is BEFORE "next" cluster of 2nd t.c.
c
if(layer_dir_2*(layer_2_n-layer_1_c) .ge. 0) then
c
c Layer condition is satisfied. 1st t.c. is in between 2
c adjacent clusters of 2nd t.c. Try to insert it between them.
c
insert_status = insert_inbetween(
& layer_2_c,layer_2_n,layer_1_c,
& inddat_dprs_2,nhits2,inddat_dprs_1,nhits1,
& hitn_2_c_f,hitn_2_c_l,hitn_2_n_f,hitn_2_n_l,
& hitn_1_c_f,hitn_1_c_l,
& num_w_saved,same_layer_f)
if(insert_status.lt.0) then
c
c Insertion failed. Only if 2nd t.c. changes its movement direction it
c can have another pair of adjacent clusters sutisfying the above criteria
c
if(iw(inddat_dprs_2-DPRCEN+DPRLRG) .gt. 0) go to 111
else
c
c Insertion scceeded. Enlarge DPRS of the 2nd t.c., modify it and exit
c
if(same_layer_f) then
c
c "next" cluster of 2nd t.c. has been saved to merging buffer
c
call modify_DPRS_DHRE(
& itr2,nhits2,inddat_dprs_2,hitn_2_c_f-1,
& itr2,nhits2,inddat_dprs_2,hitn_2_n_l+1,
& num_w_saved,itr1,nhits1,inddat_dprs_1)
else
c
c "next" cluster of 2nd t.c. has NOT been saved to merging buffer
c
call modify_DPRS_DHRE(
& itr2,nhits2,inddat_dprs_2,hitn_2_c_f-1,
& itr2,nhits2,inddat_dprs_2,hitn_2_n_f,
& num_w_saved,itr1,nhits1,inddat_dprs_1)
endif
go to 111
endif
else
c
c Layer condition is not satisfied. Go to next cluster of 2nd t.c.
c Which means
c 1. redefine "next" cluster of the 2nd t.c. as "current" one
c 2. take new cluster of the 2nd t.c. as "next"
c 3. continue main loop
c
hitn_2_c_f = hitn_2_n_f
hitn_2_c_l = hitn_2_n_l
layer_2_c = layer_2_n
hitn_2_n_f = hitn_2_c_l+1
if(hitn_2_n_f.le.nhits2) then
hitn_2_n_l = get_layer_cluster(
& nhits2,hitn_2_n_f,
& inddat_dprs_2,layer_2_n)
endif
endif
enddo
else if(layer_dir_2*(layer_2_c-layer_1_c) .eq. 0) then
c
c Both t.c. have clusters in the same layer. Combine this clusters in one
c
n_comb_wires=combine_clusters(
& hitn_1_c_f,hitn_1_c_l,inddat_dprs_1,nhits1,
& hitn_2_c_f,hitn_2_c_l,inddat_dprs_2,nhits2,
& combine_buffer)
if(n_comb_wires.gt.0) then
c
c Combining succeeded.
c Save combined cluster into merging buffer
c
cycle=1
do while(cycle.le.n_comb_wires)
merge_buff(cycle) = combine_buffer(cycle)
cycle = cycle + 1
enddo
c
c Enlarge DPRS of the 2nd t.c., modify it and exit
c
call modify_DPRS_DHRE(
& itr2,nhits2,inddat_dprs_2,0,
& itr2,nhits2,inddat_dprs_2,hitn_2_n_f,
& num_w_saved,itr1,nhits1,inddat_dprs_1)
go to 111
else
c
c Combining failed. Only if 2nd t.c. changes its movement direction it
c can have another pair of adjacent clusters sutisfying the above criteria
c
if(iw(inddat_dprs_2-DPRCEN+DPRLRG) .gt. 0) go to 111
endif
c else
c
c 1st t.c. is BEFORE "current" cluster of 2nd t.c., but within LAYER_GAP
c (otherwise this routine is never called). Check azimuth direction and gap
c and either merge these t.c. or leave this routine.
c It's not an overlaping of t.c., so at best we can have only one connection
c and only one z-measurement. It's not sufficient to start track fit and in
c this sence the merged t.c. will not be really 3-d.
c For a moment such t.c. are net merged
c
endif
111 return
end
C===========================================================================
subroutine merge_2_tc(
& itr1,inddat_dprs_1i,itr2,inddat_dprs_2i)
$$IMPLICIT
integer itr1,inddat_dprs_1i,itr2,inddat_dprs_2i
C-----------------------------------------------------------------------
C-
C- Created 20-May-1996 Alexander Andryakov
C-
C- Purpose: Merge 2 t.c.
C-
C- Input :
C- itr1 \ DPRS bank numbers of both t.c.
C- itr2 /
C- inddat_dprs_1i \ pointer to the beginning of the data block of
C- inddat_dprs_2i / DPRS banks of 1st and 2nd t.c. respectevly
C- Output :
C----------------------------------------------------------------------
$$INCLUDE 'K$INC:BCS.INC'
$$INCLUDE 'YBOS$LIBRARY:ERRCOD.INC'
$$INCLUDE 'K$ITRK:DPRS.INC'
$$INCLUDE 'K$ITRK:DHRE.INC'
$$INCLUDE 'K$ITRK:DCMERGINT.INC'
c
INTEGER get_layer_cluster,save_clust2buf,combine_clusters,
& insert_inbetween, modify_DPRS_DHRE
LOGICAL check_azim_gap
c
INTEGER
& nhits1,nhits2,view_n1,view_n2,
& inddat_dprs_1,inddat_dprs_2,
& hitn_1_c_f,hitn_1_c_l,hitn_1_n_f,hitn_1_n_l,
& hitn_2_c_f,hitn_2_c_l,hitn_2_n_f,hitn_2_n_l,
& layer_1_c,layer_1_n,layer_2_c,layer_2_n,
& num_w_saved,insert_status,layer_dir_1,layer_dir_2,
& work_l,
& itr_save,nhits_save,inddat_dprs_save,hitn_save_l,
& itr_rest,nhits_rest,inddat_dprs_rest,hitn_rest_f,
& itr_drop,nhits_drop,inddat_dprs_drop,
& n_comb_wires,combine_buffer(BUFF_LENGTH),cycle
LOGICAL same_layer_f
c
c Some preparations
c
nhits1 = iw(inddat_dprs_1i+DPRPOS) + iw(inddat_dprs_1i+DPRNEG)
nhits2 = iw(inddat_dprs_2i+DPRPOS) + iw(inddat_dprs_2i+DPRNEG)
view_n1= iw(inddat_dprs_1i+DPRPOK)
view_n2= iw(inddat_dprs_2i+DPRPOK)
inddat_dprs_1 = inddat_dprs_1i + DPRCEN
inddat_dprs_2 = inddat_dprs_2i + DPRCEN
if(view_n1 .eq. 3) then
itr_rest = itr1
nhits_rest = nhits1
inddat_dprs_rest = inddat_dprs_1
else if(view_n2 .eq. 3) then
itr_rest = itr2
nhits_rest = nhits2
inddat_dprs_rest = inddat_dprs_2
endif
if(nhits1 .ge. nhits2) then
itr_save=itr1
nhits_save=nhits1
inddat_dprs_save=inddat_dprs_1
itr_drop=itr2
nhits_drop=nhits2
inddat_dprs_drop=inddat_dprs_2
else
itr_save=itr2
nhits_save=nhits2
inddat_dprs_save=inddat_dprs_2
itr_drop=itr1
nhits_drop=nhits1
inddat_dprs_drop=inddat_dprs_1
endif
hitn_save_l=0
num_w_saved=0
n_comb_wires=0
insert_status=-10
c
c Get two 1st clusters of both t.c. and define layer direction of movement
c
hitn_1_c_f = 1
hitn_1_c_l = get_layer_cluster(nhits1,hitn_1_c_f,
& inddat_dprs_1,layer_1_c)
hitn_2_c_f = 1
hitn_2_c_l = get_layer_cluster(nhits2,hitn_2_c_f,
& inddat_dprs_2,layer_2_c)
hitn_1_n_f = hitn_1_c_l+ 1
if(hitn_1_n_f .le. nhits1) then
hitn_1_n_l = get_layer_cluster(nhits1,hitn_1_n_f,
& inddat_dprs_1,layer_1_n)
layer_dir_1 = layer_1_n - layer_1_c
else
layer_dir_1 = 0
endif
hitn_2_n_f = hitn_2_c_l + 1
if(hitn_2_n_f .le. nhits2) then
hitn_2_n_l = get_layer_cluster(nhits2,hitn_2_n_f,
& inddat_dprs_2,layer_2_n)
layer_dir_2 = layer_2_n - layer_2_c
else
layer_dir_2 = 0
endif
c
c Check if both t.c. have the same layer direction (or at least their pieces)
c
if(layer_dir_2*layer_dir_1 .lt. 0) then
c
c Layer directions at the beginning of t.c. are different
c Check if one or both t.c. change layer direction of movement
c
if(iw(inddat_dprs_1i+DPRLRG) .gt. 0 .and.
& iw(inddat_dprs_2i+DPRLRG) .gt. 0) then
c
c Neither of 2 t.c. changes layer direction. T.c. cannot be merged.
c
go to 112
else
c
c At least one of t.c. changes direction. So there is a hope that its piece
c of other direction can be merge with other t.c.
c
call adjust_layer_directions(
& nhits1,inddat_dprs_1,
& hitn_1_c_f,hitn_1_c_l,layer_1_c,layer_dir_1,
& nhits2,inddat_dprs_2,
& hitn_2_c_f,hitn_2_c_l,layer_2_c,layer_dir_2)
endif
else if(layer_dir_2*layer_dir_1 .eq. 0) then
c
c At least one t.c. is 1 layer t.c..
c We do not need to go through the main merging loop. Simplear merging will
c be performed.
c Find which one and merge it into another
c
if(layer_dir_1 .eq. 0 .and. layer_dir_2 .eq. 0) then
c
c Both t.c. are single layer t.c.
c If it's not a turning region of the track, these t.c. will not produce
c a really 3-d t.c., since there can be only one connection between
c layers of different stereo views and therefore only 1 z-measurement.
c So it's pointless to merge them, apart from the cases when there are
c many single layer t.c. of one track and merging all of them can produce
c really 3-d t.c. I think it's very rear case and probably it should be
c done during one view merging (or one view t.c. extension)
c Second case when an attempt to merge 2 single layer t.c. should be done
c is the case of turning region. In this case one of t.c. can have a
c big azinuth gap and represent two pieces of track. So second
c t.c. should be inserted between two, let's say azimuth, clusters of
c first one.
c I do not realize this for a moment, before check how important it is
c
go to 112
else if(layer_dir_1 .eq. 0) then
c
c 1st t.c. is a single layer t.c.
c Check if it can be merged into 2nd t.c.
c
if(layer_dir_2*(layer_2_c-layer_1_c) .le. LAYER_GAP)
c
c 1st t.c. is before the begining of 2nd t.c. but within layer gap.
c Try to merge
c
& call merge_1_into_2(
& itr1,nhits1,inddat_dprs_1,
& hitn_1_c_f,hitn_1_c_l,layer_1_c,
& itr2,nhits2,inddat_dprs_2,
& hitn_2_c_f,hitn_2_c_l,layer_2_c,
& hitn_2_n_f,hitn_2_n_l,layer_2_n)
c
c Then no matter what is the result exit
c
go to 112
else if(layer_dir_2 .eq. 0) then
c
c 2nd t.c. is a single layer t.c.
c Check if it can be merged into 2nd t.c.
c
if(layer_dir_1*(layer_1_c-layer_2_c) .le. LAYER_GAP)
c
c 2nd t.c. is before the begining of 1st t.c. but within layer gap.
c Try to merge
c
& call merge_1_into_2(
& itr2,nhits2,inddat_dprs_2,
& hitn_2_c_f,hitn_2_c_l,layer_2_c,
& itr1,nhits1,inddat_dprs_1,
& hitn_1_c_f,hitn_1_c_l,layer_1_c,
& hitn_1_n_f,hitn_1_n_l,layer_1_n)
c
c Then no matter what is the result exit
c
go to 112
endif
endif
c
c 2. Loop over wires/clusters of 1st & 2nd t.c. (hit lists merging)
c
do while (hitn_1_c_f .le. nhits1 .and.
& hitn_2_c_f .le. nhits2)
c
c The t.c. with > 1 layer are used in this loop
c It's also supposed that they are positioned at beginning of pieces of the
c same layer direction
c
c Find which of the 2 clusters of 2 t.c. has smallest R (biggest layer #)
c
if(layer_dir_1*(layer_1_c-layer_2_c) .lt. 0) then
c
c "Current" cluster of 1st t.c. is inner than that of 2nd t.c.
c Take next layer cluster of 1st t.c. and check if current cluster
c of 2nd t.c. is in between.
c
hitn_1_n_f = hitn_1_c_l+1
if(hitn_1_n_f.le.nhits1) then
hitn_1_n_l = get_layer_cluster(nhits1,hitn_1_n_f,
& inddat_dprs_1,layer_1_n)
layer_dir_1=layer_1_n-layer_1_c
else
c
c One layer cluster remains in 1st t.c. We have almost done the job.
c Check what happens with previous clusters
c
if(insert_status .eq. 0) then
c
c Insertion/Combining of most previous clusters succeeded.
c Merging is in progress. Check layer and azimuth gaps between two current
c clusters
c
if(abs(layer_1_c - layer_2_c) .le. LAYER_GAP) then
c
c Layer gap between "current" clusters is ok
c Check now azimuthal distance between last wire of 1st t.c. cluster
c and first wire of 2nd t.c. cluster Azimuthal gap is not constant
c and depends on the layer numbers.
c
if(check_azim_gap(hitn_1_c_l,inddat_dprs_1,
& hitn_2_c_f,inddat_dprs_2)) then
c
c Azimutal gap is ok. Connection is established.
c We should save into merging buffer "current" cluster of 1st t.c and copy
c into modified DPRS the rest of 2nd t.c. starting with hitn_2_c_f
c
num_w_saved=save_clust2buf(
& inddat_dprs_1,hitn_1_c_f,hitn_1_c_l,
& num_w_saved)
hitn_rest_f = hitn_2_c_f
itr_rest = itr2
nhits_rest = nhits2
inddat_dprs_rest = inddat_dprs_2
else
c
c Azimutal gap requirement is NOT satisfied.
c Connection is NOT established.
c The "rest" for modifying DPRS is "current" cluster of 1st t.c.
c
hitn_rest_f = hitn_1_c_f
itr_rest = itr1
nhits_rest = nhits1
inddat_dprs_rest = inddat_dprs_1
endif
else
c
c Layer gap requirement is NOT satisfied.
c Connection is NOT established
c The "rest" for modifying DPRS is "current" cluster of 1st t.c.
c
hitn_rest_f = hitn_1_c_f
itr_rest = itr1
nhits_rest = nhits1
inddat_dprs_rest = inddat_dprs_1
endif
c else
c
c The merging failed shortly before these two clusters.
c In the present code it's impossible, since this case is treated at the
c moment when insertiona fails.
c
c if(num_w_saved.gt.0) then
c
c Some pieces of t.c. have been merged. Update the DPRS end exit
c
c endif
endif
c
c Modify DPRS and exit
c
go to 111
endif
c
c We have "next" cluster of 1st t.c.. Check it for second connection
c
if(layer_dir_1*(layer_1_n-layer_2_c) .gt. 0) then
c
c Layer condition is satisfied. Cluster of 2nd t.c. is in between 2
c adjacent clusters of 1st t.c. Try to insert it between them.
c
insert_status = insert_inbetween(
& layer_1_c,layer_1_n,layer_2_c,
& inddat_dprs_1,nhits1,inddat_dprs_2,nhits2,
& hitn_1_c_f,hitn_1_c_l,hitn_1_n_f,hitn_1_n_l,
& hitn_2_c_f,hitn_2_c_l,
& num_w_saved,same_layer_f)
if(insert_status .lt. 0 .and. num_w_saved .gt. 0) then
c
c Insertion failed and something has been already saved into merging
c buffer this means that t.c. are merged up to this point. Stop merging.
c Modify DPRS by copying merging buffer and the rest of t.c., "current"
c cluster of which was used as "next" at previous insertion step (in
c other words this cluster has been already checked on having connection)
c The "rest" quantaies have been already defined at the previous merging
c stage
c
go to 111
endif
c
c Nothing has been already saved into merging buffer or insertion succeeded
c continue procedure.
c
if(same_layer_f) then
c
c The only important thing after insertion attempt is whether
c the "next" cluster of 1st t.c. and the last one of possible
c supercluster of the 2nd t.c. are in the same layer or not.
c
c The "next" cluster of 1st t.c. and the last one of possible
c supercluster of the 2nd t.c. are IN THE SAME layer,
c
c 1. take new clusters of both t.c. as "current"
c
hitn_1_c_f = hitn_1_n_l+1
if(hitn_1_c_f.le.nhits1) then
hitn_1_c_l = get_layer_cluster(
& nhits1,hitn_1_c_f,
& inddat_dprs_1,work_l)
layer_dir_1=work_l-layer_1_c
layer_1_c=work_l
endif
hitn_2_c_f = hitn_2_c_l+1
if(hitn_2_c_f.le.nhits2) then
hitn_2_c_l = get_layer_cluster(
& nhits2,hitn_2_c_f,
& inddat_dprs_2,work_l)
layer_dir_2=work_l-layer_2_c
layer_2_c=work_l
endif
c
c There is a question in this case : how to proceed with modification of
c DPRS banks if merging stops at this step (attempt to merge next clusters
c will fail). Namely,the hits of which of two t.c. should be copied as the
c rest into the DPRS bank of new 3d t.c. In contrast to real insertion of
c the wire cluster between to others, in combining the clusters are in the
c same layer, so there is no preference (in case of combining failure)
c which of 2 t.c. shoud be taken to be the end of new 3d candidate.
c I would like to stress that the logic of this merging scheme imply that
c one of t.c. is a 3d t.c. and another is 1-view t.c.. It seems to me more
c reasonable to keep alive already merged 3d and kill the rest of 1-view
c t.c.
c Of course, this considerations have a sence only if both t.c.
c have "next" clusters. Otherwise the "rest" t.c. is obvious.
c
if(hitn_1_c_f.le.nhits1 .and. hitn_2_c_f.le.nhits2) then
c
c Both t.c. have NOT reacehed their ends
c
if(view_n1 .eq. 3) then
c
c 1st t.c. is 3d one. The "rest" for modifying DPRS starts at
c "current" cluster of 1st t.c.
c
hitn_rest_f = hitn_1_c_f
itr_rest = itr1
nhits_rest = nhits1
inddat_dprs_rest = inddat_dprs_1
else if(view_n2 .eq. 3) then
c
c 2nd t.c. is 3d one. The "rest" for modifying DPRS starts at
c "current" cluster of 2nd t.c.
c
hitn_rest_f = hitn_2_c_f
itr_rest = itr2
nhits_rest = nhits2
inddat_dprs_rest = inddat_dprs_2
endif
else if(hitn_1_c_f.le.nhits1) then
c
c 2nd t.c. ends at "current" cluster. The "rest" for modifying DPRS
c starts at "current" cluster of 1st t.c.
c
hitn_rest_f = hitn_1_c_f
itr_rest = itr1
nhits_rest = nhits1
inddat_dprs_rest = inddat_dprs_1
else if(hitn_2_c_f.le.nhits2) then
c
c 1st t.c. ends at "current" cluster. The "rest" for modifying DPRS
c starts at "current" cluster of 2nd t.c.
c
hitn_rest_f = hitn_2_c_f
itr_rest = itr2
nhits_rest = nhits2
inddat_dprs_rest = inddat_dprs_2
else
c
c Both t.c. have reacehed their ends. There is NO any "rest's" to
c save in modifying DPRS. To be coinsistent with the logic of the
c modify_DPRS_DHRE routine let's use as the rest 1st t.c. and define
c hitn_rest_f to be out of the hit range of the 1st t.c.
c
itr_rest = itr1
nhits_rest = nhits1
hitn_rest_f = nhits_rest +1
inddat_dprs_rest = inddat_dprs_1
endif
else
c
c The "next" cluster of 1st t.c. and the last one of possible
c supercluster of the 2nd t.c. are NOT in the same layer,
c 1. redefine "next" cluster of the 1st t.c. as "current" one
c 2. define "current" cluster of 1st t.c. as the beginning of hit list
c to be saved after merging buffer in modified DPRS
c 3. take new cluster of the 2nd t.c. as "current"
c
hitn_1_c_f = hitn_1_n_f
hitn_1_c_l = hitn_1_n_l
layer_1_c = layer_1_n
c
hitn_rest_f = hitn_1_c_f
itr_rest = itr1
nhits_rest = nhits1
inddat_dprs_rest = inddat_dprs_1
c
hitn_2_c_f = hitn_2_c_l+1
if(hitn_2_c_f.le.nhits2) then
hitn_2_c_l = get_layer_cluster(
& nhits2,hitn_2_c_f,
& inddat_dprs_2,work_l)
layer_dir_2=work_l-layer_2_c
layer_2_c=work_l
endif
endif
else if(layer_dir_1*(layer_1_n-layer_2_c) .lt. 0) then
c
c "next" layer cluster of 1st t.c. is also before than "current" cluster
c of 2nd t.c.. Extend "current" cluster of 1st t.c. to supercluster by
c adding to it its "next" cluster and continue loop over wires/clusters
c
c This scheme should be improved for t.c. which change their movement
c direction
c
hitn_1_c_l = hitn_1_n_l
layer_1_c = layer_1_n
else if(layer_1_n .eq. layer_2_c) then
c
c "next" cluster of 1st t.c. and "current" one of 2nd t.c. are in the same
c layer.
c 1. Save into merging buffer "current" cluster of 1st t.c.. (only
c in case when isertion/combining of the previous clusters has
c succeeded)
c 2. Redefine "next" cluster of 1st t.c. as "current" and go to the
c beginning of cluster loop. (the case of combining two clusters of the
c same layer will be treated correctly).
c 3. The "rest" for modifying DPRS is "next" cluster of 1st t.c.
c
if(insert_status .ge. 0 .or. n_comb_wires .gt.0) then
num_w_saved=save_clust2buf(
& inddat_dprs_1,hitn_1_c_f,hitn_1_c_l,
& num_w_saved)
else
c
c Special case of combining. One of t.c. obviously starts closer to
c interaction region than another. We should redefine "to be saved"
c and "to be dropped" t.c. in a such way that one with greater
c layer number of the first cluster to be saved and another to be
c dropped.
c
hitn_save_l= hitn_1_c_l
itr_save = itr1
nhits_save = nhits1
inddat_dprs_save = inddat_dprs_1
itr_drop=itr2
nhits_drop=nhits2
inddat_dprs_drop=inddat_dprs_2
endif
c
hitn_1_c_f = hitn_1_n_f
hitn_1_c_l = hitn_1_n_l
layer_1_c = layer_1_n
c
hitn_rest_f = hitn_1_n_f
itr_rest = itr1
nhits_rest = nhits1
inddat_dprs_rest = inddat_dprs_1
endif
else if(layer_dir_2*(layer_2_c-layer_1_c) .lt. 0) then
c
c "Current" cluster of 2nd t.c. is inner than that of 1st t.c.
c Take next layer cluster of 2nd t.c. and check if current cluster
c of 1st t.c. is in between.
c
hitn_2_n_f = hitn_2_c_l+1
if(hitn_2_n_f.le.nhits2) then
hitn_2_n_l = get_layer_cluster(nhits2,hitn_2_n_f,
& inddat_dprs_2,layer_2_n)
layer_dir_2=layer_2_n-layer_2_c
else
c
c One layer cluster remains in 2nd t.c. We have almost done the job.
c Check what happens with previous clusters
c
if(insert_status .eq. 0) then
c
c Insertion/Combining of most previous clusters succeeded.
c Merging is in progress. Check layer and azimuth gaps between two current
c clusters
c
if(abs(layer_1_c - layer_2_c) .le. LAYER_GAP) then
c
c Layer gap between "current" clusters is ok
c Check now azimuthal distance between last wire of 2nd t.c. cluster
c and first wire of 1st t.c. cluster Azimuthal gap is not constant
c and depends on the layer numbers.
c
if(check_azim_gap(hitn_2_c_l,inddat_dprs_2,
& hitn_1_c_f,inddat_dprs_1)) then
c
c Azimutal gap is ok. Connection is established.
c We should save into merging buffer "current" cluster of 2nd t.c and copy
c into modified DPRS the rest of 1st t.c. starting with hitn_1_c_f
c
num_w_saved=save_clust2buf(
& inddat_dprs_2,hitn_2_c_f,hitn_2_c_l,
& num_w_saved)
hitn_rest_f = hitn_1_c_f
itr_rest = itr1
nhits_rest = nhits1
inddat_dprs_rest = inddat_dprs_1
else
c
c Azimutal gap requirement is NOT satisfied.
c Connection is NOT established.
c The "rest" for modifying DPRS is "current" cluster of 2nd t.c.
c
hitn_rest_f = hitn_2_c_f
itr_rest = itr2
nhits_rest = nhits2
inddat_dprs_rest = inddat_dprs_2
endif
else
c
c Layer gap requirement is NOT satisfied.
c Connection is NOT established
c The "rest" for modifying DPRS is "current" cluster of 2nd t.c.
c
hitn_rest_f = hitn_2_c_f
itr_rest = itr2
nhits_rest = nhits2
inddat_dprs_rest = inddat_dprs_2
endif
c else
c
c The merging failed shortly before these two clusters.
c In the present code it's impossible, since this case is treated at the
c moment when insertiona fails.
c
c if(num_w_saved.gt.0) then
c
c Some pieces of t.c. have been merged. Update the DPRS end exit
c
c endif
endif
c
c Modify DPRS and exit
c
go to 111
endif
c
c We have "next" cluster of 2nd t.c.. Check it for second connection
c
if(layer_dir_2*(layer_2_n-layer_1_c) .gt. 0) then
c
c Layer condition is satisfied. Cluster of 1st t.c. is in between 2
c adjacent clusters of 2nd t.c. Try to insert it between them.
c
insert_status = insert_inbetween(
& layer_2_c,layer_2_n,layer_1_c,
& inddat_dprs_2,nhits2,inddat_dprs_1,nhits1,
& hitn_2_c_f,hitn_2_c_l,hitn_2_n_f,hitn_2_n_l,
& hitn_1_c_f,hitn_1_c_l,
& num_w_saved,same_layer_f)
if(insert_status .lt. 0 .and. num_w_saved .gt. 0) then
c
c Insertion failed and something has been already saved into merging
c buffer this means that t.c. are merged up to this point. Stop merging.
c Modify DPRS by copying merging buffer and the rest of t.c., "current"
c cluster of which was used as "next" at previous insertion step (in
c other words this cluster has been already checked on having connection)
c
c The "rest" quantaies have been already defined at the previous merging
c stage
c
go to 111
endif
if(same_layer_f) then
c
c The only important thing after insertion attempt is whether
c the "next" cluster of 2nd t.c. and the last one of possible
c supercluster of the 1st t.c. are in the same layer or not.
c
c The "next" cluster of 2nd t.c. and the last one of possible
c supercluster of the 1st t.c. are IN THE SAME layer,
c 1. take new clusters of both t.c. as "current"
c
hitn_2_c_f = hitn_2_n_l+1
if(hitn_2_c_f.le.nhits2) then
hitn_2_c_l = get_layer_cluster(
& nhits2,hitn_2_c_f,
& inddat_dprs_2,work_l)
layer_dir_2=work_l-layer_2_c
layer_2_c=work_l
endif
hitn_1_c_f = hitn_1_c_l+1
if(hitn_1_c_f.le.nhits1) then
hitn_1_c_l = get_layer_cluster(
& nhits1,hitn_1_c_f,
& inddat_dprs_1,work_l)
layer_dir_1=work_l-layer_1_c
layer_1_c=work_l
endif
c
if(hitn_1_c_f.le.nhits1 .and. hitn_2_c_f.le.nhits2) then
c
c Both t.c. have NOT reacehed their ends
c
if(view_n1 .eq. 3) then
c
c 1st t.c. is 3d one. The "rest" for modifying DPRS starts at
c "current" cluster of 1st t.c.
c
hitn_rest_f = hitn_1_c_f
itr_rest = itr1
nhits_rest = nhits1
inddat_dprs_rest = inddat_dprs_1
else if(view_n2 .eq. 3) then
c
c 2nd t.c. is 3d one. The "rest" for modifying DPRS starts at
c "current" cluster of 2nd t.c.
c
hitn_rest_f = hitn_2_c_f
itr_rest = itr2
nhits_rest = nhits2
inddat_dprs_rest = inddat_dprs_2
endif
else if(hitn_1_c_f.le.nhits1) then
c
c 2nd t.c. ends at "current" cluster. The "rest" for modifying DPRS
c starts at "current" cluster of 1st t.c.
c
hitn_rest_f = hitn_1_c_f
itr_rest = itr1
nhits_rest = nhits1
inddat_dprs_rest = inddat_dprs_1
else if(hitn_2_c_f.le.nhits2) then
c
c 1st t.c. ends at "current" cluster. The "rest" for modifying DPRS
c starts at "current" cluster of 2nd t.c.
c
hitn_rest_f = hitn_2_c_f
itr_rest = itr2
nhits_rest = nhits2
inddat_dprs_rest = inddat_dprs_2
else
c
c Both t.c. have reacehed their ends. There is NO any "rest's" to
c save in modifying DPRS. To be coinsistent with the logic of the
c modify_DPRS_DHRE routine let's use as the rest 1st t.c. and define
c hitn_rest_f to be out of the hit range of the 1st t.c.
c
itr_rest = itr1
nhits_rest = nhits1
hitn_rest_f = nhits_rest +1
inddat_dprs_rest = inddat_dprs_1
endif
else
c
c The "next" cluster of 2nd t.c. and the last one of possible
c supercluster of the 1st t.c. are NOT in the same layer,
c 1. redefine "next" cluster of the 2nd t.c. as "current" one
c 2. define "current" cluster of 2nd t.c. as the beginning of hit list
c to be saved after merging buffer in modified DPRS
c 3. take new cluster of the 1st t.c. as "current"
c
hitn_2_c_f = hitn_2_n_f
hitn_2_c_l = hitn_2_n_l
layer_2_c = layer_2_n
c
hitn_rest_f = hitn_2_c_f
itr_rest = itr2
nhits_rest = nhits2
inddat_dprs_rest = inddat_dprs_2
c
hitn_1_c_f = hitn_1_c_l+1
if(hitn_1_c_f.le.nhits1) then
hitn_1_c_l = get_layer_cluster(
& nhits1,hitn_1_c_f,
& inddat_dprs_1,work_l)
layer_dir_1=work_l-layer_1_c
layer_1_c=work_l
endif
endif
else if(layer_dir_2*(layer_2_n-layer_1_c) .lt. 0) then
c
c "next" layer cluster of 2nd t.c. is also before than "current" cluster
c of 1st t.c.. Extend "current" cluster of 2nd t.c. to supercluster by
c adding to it its "next" cluster and continue loop over wires/clusters
c
c This scheme should be improved for t.c. which change their movement
c direction
c
hitn_2_c_l = hitn_2_n_l
layer_2_c = layer_2_n
else if(layer_2_n .eq. layer_1_c) then
c
c "next" cluster of 2nd t.c. and "current" one of 1st t.c. are in the same
c layer.
c 1. Save into merging buffer "current" cluster of 2nd t.c..(only
c in case when isertion/combining of the previous clusters has
c succeeded)
c 2. Redefine "next" cluster of 2nd t.c. as "current" and go to the
c beginning of cluster loop. (the case of combining two clusters of the
c same layer will be treated correctly).
c 3. The "rest" for modifying DPRS is "next" cluster of 2nd t.c.
c
if(insert_status .ge. 0 .or. n_comb_wires .gt.0) then
num_w_saved=save_clust2buf(
& inddat_dprs_2,hitn_2_c_f,hitn_2_c_l,
& num_w_saved)
else
c
c Special case of combining. One of t.c. obviously starts closer to
c interaction region than another. We should redefine "to be saved"
c and "to be dropped" t.c. in a such way that one with greater
c layer number of the first cluster to be saved and another to be
c dropped.
c
hitn_save_l= hitn_2_c_l
itr_save = itr2
nhits_save = nhits2
inddat_dprs_save = inddat_dprs_2
itr_drop=itr1
nhits_drop=nhits1
inddat_dprs_drop=inddat_dprs_1
endif
c
hitn_2_c_f = hitn_2_n_f
hitn_2_c_l = hitn_2_n_l
layer_2_c = layer_2_n
c
hitn_rest_f = hitn_2_n_f
itr_rest = itr2
nhits_rest = nhits2
inddat_dprs_rest = inddat_dprs_2
endif
else
c
c "Current" clusters of both t.c. are in the same layer. They should be
c combined.
c First, define azimuth direction, on the base of layer direction and sign
c of curvature
c
if(layer_dir_1 .ne. 0) then
phi_dir_cur_1=sign(1.,layer_dir_1*rw(inddat_dprs_1i+DPRCUR))
else
phi_dir_cur_1=0
endif
if(layer_dir_2 .ne. 0) then
phi_dir_cur_2=sign(1.,layer_dir_2*rw(inddat_dprs_2i+DPRCUR))
else
phi_dir_cur_2=0
endif
c
c Perform combining
c
n_comb_wires=combine_clusters(
& hitn_1_c_f,hitn_1_c_l,inddat_dprs_1,nhits1,
& hitn_2_c_f,hitn_2_c_l,inddat_dprs_2,nhits2,
& combine_buffer)
if(n_comb_wires .lt. 0 .and. num_w_saved .gt. 0)
c
c Combining failed and something has been already saved into merging
c buffer this means that t.c. are merged up to this point. Stop merging.
c Modify DPRS by copying merging buffer and the rest of t.c., "current"
c cluster of which was used as "next" at previous insertion step (in
c other words this cluster has been already checked on having connection)
c The "rest" quantaies have been already defined at the previous merging
c stage
c
& go to 111
c
c The clusters of the same layer have been succsessfully combined.
c Copy the combining buffer into merging one
c
cycle=1
do while(cycle.le.n_comb_wires .and.
& num_w_saved .lt. BUFF_LENGTH)
num_w_saved = num_w_saved +1
merge_buff(num_w_saved) = combine_buffer(cycle)
cycle = cycle + 1
enddo
c
c Take new clusters of both t.c. as "current"
c
hitn_1_c_f = hitn_1_c_l+1
if(hitn_1_c_f.le.nhits1) then
hitn_1_c_l = get_layer_cluster(
& nhits1,hitn_1_c_f,
& inddat_dprs_1,work_l)
layer_dir_1=work_l-layer_1_c
layer_1_c=work_l
endif
hitn_2_c_f = hitn_2_c_l+1
if(hitn_2_c_f.le.nhits2) then
hitn_2_c_l = get_layer_cluster(
& nhits2,hitn_2_c_f,
& inddat_dprs_2,work_l)
layer_dir_2=work_l-layer_2_c
layer_2_c=work_l
endif
if(hitn_1_c_f.le.nhits1 .and. hitn_2_c_f.le.nhits2) then
c
c Both t.c. have NOT reacehed their ends
c
if(view_n1 .eq. 3) then
c
c 1st t.c. is 3d one. The "rest" for modifying DPRS starts at
c "current" cluster of 1st t.c.
c
hitn_rest_f = hitn_1_c_f
itr_rest = itr1
nhits_rest = nhits1
inddat_dprs_rest = inddat_dprs_1
else if(view_n2 .eq. 3) then
c
c 2nd t.c. is 3d one. The "rest" for modifying DPRS starts at
c "current" cluster of 2nd t.c.
c
hitn_rest_f = hitn_2_c_f
itr_rest = itr2
nhits_rest = nhits2
inddat_dprs_rest = inddat_dprs_2
endif
else if(hitn_1_c_f.le.nhits1) then
c
c 2nd t.c. ends at "current" cluster. The "rest" for modifying DPRS
c starts at "current" cluster of 1st t.c.
c
hitn_rest_f = hitn_1_c_f
itr_rest = itr1
nhits_rest = nhits1
inddat_dprs_rest = inddat_dprs_1
else if(hitn_2_c_f.le.nhits2) then
c
c 1st t.c. ends at "current" cluster. The "rest" for modifying DPRS
c starts at "current" cluster of 2nd t.c.
c
hitn_rest_f = hitn_2_c_f
itr_rest = itr2
nhits_rest = nhits2
inddat_dprs_rest = inddat_dprs_2
else
c
c Both t.c. have reacehed their ends. There is NO any "rest's" to
c save in modifying DPRS. To be coinsistent with the logic of the
c modify_DPRS_DHRE routine let's use as the rest 1st t.c. and define
c hitn_rest_f to be out of the hit range of the 1st t.c.
c
itr_rest = itr1
nhits_rest = nhits1
hitn_rest_f = nhits_rest +1
inddat_dprs_rest = inddat_dprs_1
endif
endif
enddo
c
c Modify DPRS on conditions established during the merging procedure
c
111 call modify_DPRS_DHRE(
& itr_save,nhits_save,inddat_dprs_save,hitn_save_l,
& itr_rest,nhits_rest,inddat_dprs_rest,hitn_rest_f,
& num_w_saved,itr_drop,nhits_drop,inddat_dprs_drop)
112 return
end
C===========================================================================
subroutine modify_DPRS_DHRE(
& itr1,nhits1,inddat_dprs_1,hitn_1_l,
& itr2,nhits2,inddat_dprs_2,hitn_2_f,
& num_w_saved,itr_drop,nhits_drop,inddat_dprs_drop)
$$IMPLICIT
integer
& itr1,nhits1,inddat_dprs_1,hitn_1_l,
& itr2,nhits2,inddat_dprs_2,hitn_2_f,
& num_w_saved,itr_drop,nhits_drop,inddat_dprs_drop
C-----------------------------------------------------------------------
C-
C- Created 4-Jun-1996 Alexander Andryakov
C-
C- Purpose: Modify DPRS bank of itr1 t.c. to include merging results.
C- It's supposed that itr1 t.c. has the first hit in the new merged
C- t.c. and the itr2 t.c. contains "tail" hits. In particular case
C- they can be the same t.c.
C- Modify DHRE bank to put 0 into DHRTRN word of the hits not
C- belonging to new merged chain
C- Input :
C- itr1 - DPRS bank nubmer of t.c. containing the 1st hit
C- itr2 - DPRS bank nubmer of t.c. containing "tail" hits
C- inddat_dprs_1 \ pointer to the beginning of the hit list of DPRS
C- inddat_dprs_2 / banks of both t.c.
C- nhits1 \ number of hits of both t.c.
C- nhits2 /
C- hitn_1_l - DHRE number of the last hit before merging region
C- hitn_2_f - DHRE number of the first hit after merging region
C- num_w_saved - number of hits in merging buffer
C- itr_drop - DPRS bank nubmer of t.c. to be dropped
C- inddat_dprs_drop - pointer to the beginning of the hit list of DPRS
C- bank to be dropped
C- nhits_drop - number of hits of this t.c.
C----------------------------------------------------------------------
$$INCLUDE 'K$INC:BCS.INC'
$$INCLUDE 'YBOS$LIBRARY:ERRCOD.INC'
$$INCLUDE 'K$INC:ERLEVL.INC'
$$INCLUDE 'K$ITRK:CTIDRI.INC'
c$$INCLUDE 'K$ITRK:CDTHIT.INC'
$$INCLUDE 'K$ITRK:DPRS.INC'
$$INCLUDE 'K$ITRK:DHRE.INC'
$$INCLUDE 'K$ITRK:DCMERGINT.INC'
C
CHARACTER*36 BNKTYP
DATA BNKTYP/'2I4,9R4,2I4,R4,2I4,2R4,I4,XXX(I4,R4)'/
C
c
INTEGER BTYMAK,BCOPY,MDROP
c
INTEGER
& nhits_new,
& status,work_lng,work_ind,work_inddat,work_itr,
& new_ind,new_inddat,
& cycle,hit_cycle,layer,nneg,npos,
& inddat_dprs1_n,dhre_w_n,dhre_pos
C
LOGICAL ch_dir,first
INTEGER layer_max,layer_min,layer_dir,layer_pr,temp
c
c For a moment
c
if(num_w_saved.le.0) return
c
work_itr = 1
c
c Stupid check
c
if(itr1.eq.itr_drop) then
write (6,'('' modify_DPRS_DHRE .... itr1,itr_drop:'',2i10)')
& itr1,itr_drop
stop
endif
c
c Change DHRE to free all hits belonging to "host" t.c.
c
cycle=1
do while(cycle .le. nhits1)
dhre_w_n = iw(inddat_dprs_1 + DPRLWI*(cycle-1)+DPRIWI)
iw(inddat_dhre + dhre_e_size*(dhre_w_n-1) + DHRTRN)=0
cycle=cycle+1
enddo
c
c Change DHRE to free all hits belonging to "gifter" t.c.
c
cycle=1
do while(cycle .le. nhits_drop)
dhre_w_n = iw(inddat_dprs_drop + DPRLWI*(cycle-1)+DPRIWI)
iw(inddat_dhre + dhre_e_size*(dhre_w_n-1) + DHRTRN)=0
cycle=cycle+1
enddo
c
c Copy into working DPRS bank the DPRS bank of the t.c. that contains first
c hit
c
status=BCOPY(iw,'DPRS',itr1,'WORK',work_itr,work_ind,work_inddat)
if(status .ne. YESUCC) then
c
c DPRS bank of 1st t.c. cannot be copied to working DPRS
c This event should be reprocessed
c
call ERLOGR('modify_DPRS_DHRE',ERSEVR,0,0,
+ 'Bank DPRS cannot be copied')
return
endif
c
c Calculate new number of hits
c
nhits_new=hitn_1_l +
& num_w_saved +
& nhits2-hitn_2_f+1
if(nhits_new.gt.999) nhits_new=999
write (BNKTYP(27:29),'(I3)') nhits_new
c
c Enlarge working DPRS bank
c
status = BTYMAK(iw,'WORK',work_itr,BNKTYP,
& work_lng,work_ind,work_inddat)
if(status.ne.YESUCC) then
c
c Working DPRS bank cannot be enlarged due to lack of space and has
c been lost. This event should be reprocessed
c
call ERLOGR('modify_DPRS_DHRE',ERSEVR,0,0,
+ 'Bank WORK cannot be enlarged')
return
endif
c Initialisation for t.c. area definition
c
layer_max=-1
layer_min=1000
ch_dir=.FALSE.
first=.TRUE.
layer_dir=0
c
c First, count the number of + and - wires before the merged part
c and assign new t.c. number to these hits in DHRE
c
work_inddat = work_inddat + iw(work_inddat+DPRHDS)
inddat_dprs1_n = work_inddat + DPRCEN
nneg=0
npos=0
hit_cycle=1
do while(hit_cycle .le. hitn_1_l .and.
& hit_cycle .le. nhits_new)
dhre_w_n = iw(inddat_dprs1_n+DPRLWI*(hit_cycle-1)+DPRIWI)
dhre_pos = inddat_dhre + dhre_e_size*(dhre_w_n-1)
iw(dhre_pos + DHRTRN)=itr1
layer = iw(dhre_pos + DHRSLR)
if(TIDST(layer).gt.0) then
npos=npos+1
else if(TIDST(layer).lt.0) then
nneg=nneg+1
endif
hit_cycle = hit_cycle + 1
c
c Redefine layer boundary of the piece of t.c.
c
if(.not.first) then
if(layer_dir.ne.0 .and. .not.ch_dir)
& ch_dir = ch_dir .or.
& layer_dir*(layer - layer_pr) .lt. 0
layer_dir = layer - layer_pr
else
first=.FALSE.
endif
layer_pr=layer
if(layer .lt. layer_min) layer_min = layer
if(layer .gt. layer_max) layer_max = layer
enddo
c
c Copy merging buffer into enlarged working DPRS starting from hit
c next to hitn_1_l and count the number of + and - wires in this sample
c
cycle=1
do while(cycle .le. num_w_saved .and.
& hit_cycle .le. nhits_new)
iw(inddat_dprs1_n+DPRLWI*(hit_cycle-1)+DPRIWI)=
& merge_buff(cycle)
dhre_pos = inddat_dhre + dhre_e_size*(merge_buff(cycle)-1)
iw(dhre_pos + DHRTRN)=itr1
layer = iw(dhre_pos + DHRSLR)
if(TIDST(layer).gt.0) then
npos=npos+1
else if(TIDST(layer).lt.0) then
nneg=nneg+1
endif
cycle=cycle+1
hit_cycle=hit_cycle+1
c
c Redefine layer boundary of the piece of t.c.
c
if(.not.first) then
if(layer_dir.ne.0 .and. .not.ch_dir)
& ch_dir = ch_dir .or.
& layer_dir*(layer - layer_pr) .lt. 0
layer_dir = layer - layer_pr
else
first=.FALSE.
endif
layer_pr=layer
if(layer .lt. layer_min) layer_min = layer
if(layer .gt. layer_max) layer_max = layer
enddo
c
c Copy the rest of DPRS bank of that t.c., which has "tail" hits
c And count the number of + and - wires in this sample
c
cycle=hitn_2_f
do while(cycle .le. nhits2 .and.
& hit_cycle .le. nhits_new)
dhre_w_n = iw(inddat_dprs_2+DPRLWI*(cycle-1)+DPRIWI)
iw(inddat_dprs1_n+DPRLWI*(hit_cycle-1)+DPRIWI)=dhre_w_n
dhre_pos = inddat_dhre + dhre_e_size*(dhre_w_n-1)
iw(dhre_pos + DHRTRN)=itr1
layer = iw(dhre_pos + DHRSLR)
if(TIDST(layer).gt.0) then
npos=npos+1
else if(TIDST(layer).lt.0) then
nneg=nneg+1
endif
cycle=cycle+1
hit_cycle=hit_cycle+1
c
c Redefine layer boundary of the piece of t.c.
c
if(.not.first) then
if(layer_dir.ne.0 .and. .not.ch_dir)
& ch_dir = ch_dir .or.
& layer_dir*(layer - layer_pr) .lt. 0
layer_dir = layer - layer_pr
else
first=.FALSE.
endif
layer_pr=layer
if(layer .lt. layer_min) layer_min = layer
if(layer .gt. layer_max) layer_max = layer
enddo
c*************************************************** debug only
c
c Stupid debug check
c
if(nneg+npos.ne.nhits_new) then
write (6,'('' merge_1_into_2 ....'',
& ''nneg, npos, nhits_new:'',3i10)')
& nneg, npos, nhits_new
stop
endif
c*************************************************** debug only
c
c Hit list has been modified.
c Modify now kinematic variables. For a moment only curvature and
c quality are modified.
c There is a small question how to modify curvature.
c In case of complete merging: nhits1 + nhits_drop = npos + nneg the
c idea is very simple curvatures are weightted according to the number
c of contributed hits.
c In case of incomplete merging: nhits1 + nhits_drop > npos + nneg
c it's not very easy to calculate the number of contributed hits. In
c this case simple avareging of 2 curvatures is performed if the
c shortest t.c. is at least 1/3 of the other one. If it's even shorter
c then the curvature of the longest t.c. is taken.
c
if(nhits1 + nhits_drop .eq. npos + nneg) then
rw(work_inddat+DPRCUR) = (nhits1*rw(work_inddat+DPRCUR) +
& nhits_drop*rw(inddat_dprs_drop-DPRCEN+DPRCUR))/
& (npos+nneg)
else
if(nhits1.gt.nhits_drop) then
if(nhits1/nhits_drop.lt.3)
& rw(work_inddat+DPRCUR) = 0.5*(rw(work_inddat+DPRCUR) +
& rw(inddat_dprs_drop-DPRCEN+DPRCUR))
else
if(nhits_drop/nhits1.lt.3) then
rw(work_inddat+DPRCUR) = 0.5*(rw(work_inddat+DPRCUR) +
& rw(inddat_dprs_drop-DPRCEN+DPRCUR))
else
rw(work_inddat+DPRCUR) = rw(inddat_dprs_drop-DPRCEN+DPRCUR)
endif
endif
endif
c
c Modyfication to the initial azimuth angle of the t.c. momentum
c Simple averaging is not good enough due to the fact that sometimes
c very short t.c. can have absolutely wrong phi0 (even of the opsite
c sign and not within 2pi).
c
if(rw(work_inddat+DPRPHI)*rw(inddat_dprs_drop-DPRCEN+DPRPHI).gt.0)
& then
rw(work_inddat+DPRPHI) = (nhits1*rw(work_inddat+DPRPHI) +
& nhits_drop*rw(inddat_dprs_drop-DPRCEN+DPRPHI))/(npos+nneg)
else if(nhits1.lt.nhits_drop) then
rw(work_inddat+DPRPHI) = rw(inddat_dprs_drop-DPRCEN+DPRPHI)
endif
c Since "quality" is "something chi**2 like" - N.D.F., it seems to
c me more natural summ up these values instead of avareging them
c
rw(work_inddat+DPRQUA) = rw(work_inddat+DPRQUA) +
& rw(inddat_dprs_drop-DPRCEN+DPRQUA)
iw(work_inddat+DPRPOS) = npos
iw(work_inddat+DPRNEG) = nneg
iw(work_inddat+DPRPOK)=3
c
c Put new boundaries of t.c.
c
temp=abs(iw(work_inddat+DPRLRG))/10000
if(ch_dir) then
iw(work_inddat+DPRLRG)=-(temp*10000 +
+ layer_max*100 + layer_min)
else
iw(work_inddat+DPRLRG)= temp*10000 +
+ layer_max*100 + layer_min
endif
rw(work_inddat+DPRPHL)=-101.
c
c Copy working DPRS into DPRS of 1st t.c. (the bank is also enlarged
c during this copy)
c
status = BCOPY(iw,'WORK',work_itr,'DPRS',itr1,
& new_ind,new_inddat)
if(status.ne.YESUCC) then
c
c Working DPRS bank cannot be copied to 1st DPRS bank
c This event should be reprocessed
c
call ERLOGR('modify_DPRS_DHRE',ERSEVR,0,0,
+ 'Bank DPRS cannot be updated')
else
inddat_dprs_array(itr1) = new_inddat + iw(new_inddat+DPRHDS)
endif
c
c Drop one of DPRS bank
c
status = MDROP(iw,'DPRS',itr_drop)
inddat_dprs_array(itr_drop)=0
c
c A little extracalculation usefull in the FIRST merging loop
c
if(itr1 .gt. 0 .and. itr1 .lt. first_itr2) then
n_itr1_checked = n_itr1_checked +1
else if(itr1 .ge. first_itr2 .and. itr1 .le. last_itr2) then
n_itr2_checked = n_itr2_checked +1
endif
if(itr_drop .gt. 0 .and. itr_drop .lt. first_itr2) then
n_itr1_checked = n_itr1_checked +1
else if(itr_drop .ge. first_itr2 .and. itr_drop .le. last_itr2)
& then
n_itr2_checked = n_itr2_checked +1
endif
return
end
C===========================================================================
integer function save_clust2buf(
& inddat_dprs,cl_f,cl_l,
& num_w_saved)
$$IMPLICIT
integer
& inddat_dprs,cl_f,cl_l,num_w_saved
C----------------------------------------------------------------------
C-
C- Created 27-May-1996 Alexander Andryakov
C-
C- Input :
C- inddat_dprs - pointer to the beginning of the data block of
C- DPRS bank of the t.c.
C- cl_f \
C- cl_l / first & last wires in the cluster of the t.c.
C- num_w_saved - number of wires already saved in this buffer
C- Return :
C- new number of wires saved in the buffer
C----------------------------------------------------------------------
$$INCLUDE 'K$INC:BCS.INC'
$$INCLUDE 'YBOS$LIBRARY:ERRCOD.INC'
$$INCLUDE 'K$ITRK:DPRS.INC'
$$INCLUDE 'K$ITRK:DCMERGINT.INC'
c
INTEGER cycle, n_w
c
cycle = 0
n_w = cl_l-cl_f
do while (num_w_saved + cycle .lt. BUFF_LENGTH .and.
& cycle .le. n_w)
merge_buff(num_w_saved + cycle +1)=
& iw(inddat_dprs + DPRLWI*(cl_f + cycle - 1)+DPRIWI)
cycle = cycle + 1
enddo
if(n_w + 1 + num_w_saved .gt. BUFF_LENGTH) then
save_clust2buf = BUFF_LENGTH + 1
else
save_clust2buf = num_w_saved + cycle
endif
return
end
C===========================================================================
subroutine z0_ctgth
$$IMPLICIT
C----------------------------------------------------------------------
C-
C- Created 3-DEC-1993 Alexander Andryakov
C- Updated 25-Jun-1993 Alexander Andryakov (for new merging code)
C-
C- Purpose: Performs 0 approximation fit of z0 and cot(theta)
C- for each 3d t.c.
C-
C----------------------------------------------------------------------
$$INCLUDE 'K$INC:BCS.INC'
$$INCLUDE 'YBOS$LIBRARY:ERRCOD.INC'
$$INCLUDE 'K$INC:ERLEVL.INC'
$$INCLUDE 'K$ITRK:CTIDRI.INC'
$$INCLUDE 'K$ITRK:CDTHIT.INC'
$$INCLUDE 'K$ITRK:CDTCUT.INC'
$$INCLUDE 'K$ITRK:CTRACK.INC'
$$INCLUDE 'K$ITRK:DPRS.INC'
$$INCLUDE 'K$ITRK:DHRE.INC'
$$INCLUDE 'K$ITRK:DCMERGINT.INC'
c
INTEGER BLOCAT,BNEXT,BDATA
c
INTEGER dhre_w_n,dhre_w_n_p,dhre_w_p,
& status,ind_dprs,inddat_dprs,inddat_dprs_bhl,
& i_tr_ass,layer_c,last_wire,
& cycle,nhif,layer_p,n_z_meas,i_point
REAL path_rphi,zr_sum,z_sum,r_sum,r2_sum,x0_1,y0_1,
& r0_1,ri0_1,dir_x_1,dir_y_1,dir_z_1,x0_2,y0_2,
& r0_2,ri0_2,dir_x_2,dir_y_2,dir_z_2,
& ddir_x,ddir_y,ddir_z,dx_0,dy_0,a1,a2,b,c,
& denominator,wire_path_1,delta_wire_path,
& x_1,y_1,z_1,wire_path_2,
& x_2,y_2,z_2,space_points(2000,3),delta_r2,delta_r,
& avcon,cot_theta,z_0
c
c Loop over track candidates found by the pattern recognition code
c Start with the first t.c. in DPRS chain (take it as current)
c
status=BLOCAT(iw,'DPRS',-1,ind_dprs,inddat_dprs)
do while (status.eq.YESUCC.and.ind_dprs.gt.0)
c
c DPRS bank of current t.c. located succesfully
c skip the header
c
inddat_dprs=inddat_dprs+IW(inddat_dprs+DPRHDS)
inddat_dprs_bhl = inddat_dprs+DPRCEN
c
c Check first if this t.c. is 3-d t.c.
c
if(iw(inddat_dprs +DPRPOK).eq.3) then
c
c Get number of hits of located t.c.
c and make some preparations for wire loop
c
nhif=iw(inddat_dprs +DPRPOS)+iw(inddat_dprs +DPRNEG)
dhre_w_n = iw(inddat_dprs_bhl + DPRIWI)
dhre_w_p = inddat_dhre + dhre_e_size*(dhre_w_n-1)
i_tr_ass= iw(dhre_w_p + DHRTRN)
layer_c = iw(dhre_w_p + DHRSLR)
layer_p=layer_c
dhre_w_n_p=dhre_w_n
path_rphi= 0.
zr_sum = 0.
z_sum = 0.
r_sum = 0.
r2_sum = 0.
n_z_meas = 0
last_wire=0
c
c Loop over wires of located t.c.
c
do cycle=2,nhif
c
c Get the wire position in DHRE bank and define the layer number
c
dhre_w_n = iw(inddat_dprs_bhl + DPRLWI*(cycle-1) +DPRIWI)
dhre_w_p = inddat_dhre + dhre_e_size*(dhre_w_n-1)
layer_c = iw(dhre_w_p + DHRSLR)
c
c Put the track number of the first hit to all others
c
iw(dhre_w_p + DHRTRN) = i_tr_ass
c
c Find point of layer changing.
c Compare layer number of the current hit with that of the previous hit
c
if(abs(layer_p-layer_c).eq.1) then
c Two adjacent hits are in ajacent layers
c What is 3-d coordinates of these two hits?
c
x0_1=DTXH_K(dhre_w_n_p) ! previous hit is 1st
y0_1=DTYH_K(dhre_w_n_p)
r0_1=x0_1**2+y0_1**2
if(r0_1.lt.1.e-10) go to 111
ri0_1=1./sqrt(r0_1)
dir_x_1=-TIDST(layer_p)*y0_1*ri0_1
dir_y_1= TIDST(layer_p)*x0_1*ri0_1
dir_z_1= TIDSC(layer_p)
x0_2=DTXH_K(dhre_w_n) ! current one is 2nd
y0_2=DTYH_K(dhre_w_n)
r0_2=x0_2**2+y0_2**2
if(r0_2.lt.1.e-10) go to 111
ri0_2=1./sqrt(r0_2)
dir_x_2=-TIDST(layer_c)*y0_2*ri0_2
dir_y_2= TIDST(layer_c)*x0_2*ri0_2
dir_z_2= TIDSC(layer_c)
ddir_x = dir_x_2-dir_x_1
ddir_y = dir_y_2-dir_y_1
ddir_z = dir_z_2-dir_z_1
dx_0=x0_2-x0_1
dy_0=y0_2-y0_1
a1=dx_0*ddir_x +dy_0*ddir_y
a2=dx_0*dir_x_2+dy_0*dir_y_2
b=ddir_x**2 + ddir_y**2 + ddir_z**2
c=ddir_x*dir_x_2 + ddir_y*dir_y_2 + ddir_z*dir_z_2
if(abs(b-c**2).lt.1.e-10) then ! system is unsolvable
write(6,'('' wires do not cross ???'')')
go to 111 ! skip hit
endif
denominator=1./(b-c**2)
wire_path_1 = -(a1 -a2*c)*denominator
delta_wire_path = -(a2*b-a1*c)*denominator
x_1=x0_1 + dir_x_1*wire_path_1 ! 3-d coordinates 1st hit
y_1=y0_1 + dir_y_1*wire_path_1
z_1= dir_z_1*wire_path_1
wire_path_2 = wire_path_1 +delta_wire_path
x_2=x0_2 + dir_x_2*wire_path_2 ! 3-d coordinates 2nd hit
y_2=y0_2 + dir_y_2*wire_path_2
z_2= dir_z_2*wire_path_2
c write(6,'('' x,y,z,l = '',3g15.4,i5)') x_1,y_1,z_1,85-layer_p
c write(6,'('' x,y,z,l = '',3g15.4,i5)') x_2,y_2,z_2,85-layer_c
space_points(n_z_meas+1,1)=x_1
space_points(n_z_meas+1,2)=y_1
space_points(n_z_meas+1,3)=z_1
space_points(n_z_meas+2,1)=x_2
space_points(n_z_meas+2,2)=y_2
space_points(n_z_meas+2,3)=z_2
n_z_meas = n_z_meas + 2
last_wire=dhre_w_n
111 continue
endif
layer_p=layer_c
dhre_w_n_p=dhre_w_n
enddo ! loop over wires of the current track
c
c Calculate cot(theta) and Z0. (l.s.f.)
c
do i_point=2,n_z_meas
delta_r2 =
& (space_points(i_point,1)-space_points(i_point-1,1))**2 +
& (space_points(i_point,2)-space_points(i_point-1,2))**2
delta_r = sqrt(delta_r2)
path_rphi= path_rphi + delta_r
zr_sum = zr_sum + space_points(i_point,3)*path_rphi
z_sum = z_sum + space_points(i_point,3)
r_sum = r_sum + path_rphi
r2_sum = r2_sum + path_rphi**2
enddo
avcon=1./(n_z_meas-1)
denominator=r2_sum-avcon*r_sum**2
if(abs(denominator).lt.1.e-10) then ! cot(theta)=infinity
cw write(6,'('' track can. '',i3,'' is near parallel to z'')')
cw & itrack
call INERR(216,1)
else
cot_theta= (zr_sum-avcon*z_sum*r_sum)/denominator
z_0 = avcon*(z_sum-cot_theta*r_sum)
c
c Put this values into DPRS
c
c Fill z0 and cot(theta).
c It should be stressed that in present version of the code z0
c corresponds to the point of first z measurement (which is not true
c initial point of the t.c. but something corresponding to the
c intersection of the two adjacent wires in two layers with OPOSITE
c sign of stereo angle)
c
rw(inddat_dprs+DPRZCO)=z_0
rw(inddat_dprs+DPRCTT)=cot_theta
c
c Fill coordinates of the last point
c Similar to what was said on z0 is applicable here for the last z.
c
rw(inddat_dprs+DPRXLS)=DTXH_K(last_wire)
rw(inddat_dprs+DPRYLS)=DTYH_K(last_wire)
rw(inddat_dprs+DPRZLS)=z_0+cot_theta*path_rphi
endif
endif ! Check if t.c. is 3d one
c
c Locate next DPRS bank
c
status=BNEXT(iw,ind_dprs,ind_dprs)
if (status.eq.YESUCC) status=BDATA(iw,ind_dprs,inddat_dprs)
enddo ! loop over DPRS banks
call vzero(DTXH_K,NUMHIT)
call vzero(DTYH_K,NUMHIT)
END
[KLOE]
[Offline Doc]
[TRK Files]
Generated with Light on Thu Apr 8 13:00:16 MET DST 1999
.
Mail comments and suggestions.
