subroutine vtx_primary(ntrk, list, error) * >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> * Calculate primary vertex and return the tracks, the vertex and * its covariance matrix. * >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> implicit none #include "seq/clinc/anlccd.inc" #include "seq/clinc/anlclev.inc" #include "prim_vertex.inc" c Externals external kget_track_index integer kget_track_index c Calling arguments integer ntrk, list(*), error c Local Variables: integer it, ik, cdtrk, i, binp, num, binphi, err integer update logical init double precision ptcut, pcut, chicut, chisq, chidf, phideg double precision beam_pos(3), beam_sig(3), Vbeam(3,3), w(10) double precision vtx(3), Vvtx(3,3), dsigma(3) data ptcut/0.200/ data pcut/0.250/ data chicut/8./ data init/.TRUE./ c >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> if(init) then init = .FALSE. call hbook1(1, 'Number of tracks in vertex', 20, 0., 20., 0.) call hbook1(10, 'chisq/df', 50, 0., 10., 0.) call hbook1(11, 'chisq/df, p = 0.0 - 0.5', 50, 0., 10., 0.) call hbook1(12, 'chisq/df, p = 0.5 - 1.0', 50, 0., 10., 0.) call hbook1(13, 'chisq/df, p = 1.0 - 1.5', 50, 0., 10., 0.) call hbook1(14, 'chisq/df, p = 1.5 - 2.0', 50, 0., 10., 0.) call hbook1(15, 'chisq/df, p = 2.0 - 2.5', 50, 0., 10., 0.) call hbook1(16, 'chisq/df, p = 2.5 - 3.0', 50, 0., 10., 0.) call hbook1(17, 'chisq/df, p = 3.0 - 3.5', 50, 0., 10., 0.) call hbook1(18, 'chisq/df, p = 3.5 - 4.0', 50, 0., 10., 0.) call hbook1(19, 'chisq/df, p = 4.0 - 4.5', 50, 0., 10., 0.) call hbook1(20, 'chisq/df, p > 4.5 ', 50, 0., 10., 0.) call hbook1(21, 'chisq/df, phi = 0 - 45', 50, 0., 10., 0.) call hbook1(22, 'chisq/df, phi = 45 - 90', 50, 0., 10., 0.) call hbook1(23, 'chisq/df, phi = 90 -135', 50, 0., 10., 0.) call hbook1(24, 'chisq/df, phi =135 -180', 50, 0., 10., 0.) call hbook1(25, 'chisq/df, phi =180 -225', 50, 0., 10., 0.) call hbook1(26, 'chisq/df, phi =225 -270', 50, 0., 10., 0.) call hbook1(27, 'chisq/df, phi =270 -315', 50, 0., 10., 0.) call hbook1(28, 'chisq/df, phi =315 -360', 50, 0., 10., 0.) call hbook1(102, 'Vertex x error (meters)', 50, 0., 250.e-6, 0.) call hbook1(103, 'Vertex y error (meters)', 50, 0., 250.e-6, 0.) call hbook1(104, 'Vertex z error (meters)', 50, 0., 250.e-6, 0.) call hbook1(110, 'beam chisq/df all', 50, 0., 10., 0.) call hbook1(111, 'beam chisq/df tracks = 1', 50, 0., 10., 0.) call hbook1(112, 'beam chisq/df tracks = 2', 50, 0., 10., 0.) call hbook1(113, 'beam chisq/df tracks = 3', 50, 0., 10., 0.) call hbook1(114, 'beam chisq/df tracks = 4', 50, 0., 10., 0.) call hbook1(115, 'beam chisq/df tracks = 5', 50, 0., 10., 0.) call hbook1(116, 'beam chisq/df tracks = 6', 50, 0., 10., 0.) call hbook1(117, 'beam chisq/df tracks = 7', 50, 0., 10., 0.) call hbook1(118, 'beam chisq/df tracks = 8', 50, 0., 10., 0.) call hbook1(119, 'beam chisq/df tracks = 9', 50, 0., 10., 0.) call hbook1(120, 'beam chisq/df tracks > 9', 50, 0., 10., 0.) call hbook1(210, 'single chidf', 50, 0., 10., 0.) call hbook1(211, 'single chidf, p = 0.0 - 0.5', 50, 0., 10., 0.) call hbook1(212, 'single chidf, p = 0.5 - 1.0', 50, 0., 10., 0.) call hbook1(213, 'single chidf, p = 1.0 - 1.5', 50, 0., 10., 0.) call hbook1(214, 'single chidf, p = 1.5 - 2.0', 50, 0., 10., 0.) call hbook1(215, 'single chidf, p = 2.0 - 2.5', 50, 0., 10., 0.) call hbook1(216, 'single chidf, p = 2.5 - 3.0', 50, 0., 10., 0.) call hbook1(217, 'single chidf, p = 3.0 - 3.5', 50, 0., 10., 0.) call hbook1(218, 'single chidf, p = 3.5 - 4.0', 50, 0., 10., 0.) call hbook1(219, 'single chidf, p = 4.0 - 4.5', 50, 0., 10., 0.) call hbook1(220, 'single chidf, p > 4.5 ', 50, 0., 10., 0.) call hbook1(221, 'single chidf, phi = 0 - 45', 50, 0., 10., 0.) call hbook1(222, 'single chidf, phi = 45 - 90', 50, 0., 10., 0.) call hbook1(223, 'single chidf, phi = 90 -135', 50, 0., 10., 0.) call hbook1(224, 'single chidf, phi =135 -180', 50, 0., 10., 0.) call hbook1(225, 'single chidf, phi =180 -225', 50, 0., 10., 0.) call hbook1(226, 'single chidf, phi =225 -270', 50, 0., 10., 0.) call hbook1(227, 'single chidf, phi =270 -315', 50, 0., 10., 0.) call hbook1(228, 'single chidf, phi =315 -360', 50, 0., 10., 0.) endif * Get beam position, beam size and beam covariance matrix call kget_beam_position(beam_pos, beam_sig) call kget_beam_position_cov(Vbeam) num_vtx = 0 do it=1,ntrk ik = list(it) call kget_track_param(ik, w) !Put track parameters in w(1-10) cdtrk = kget_track_index(ik) !CD track index if(kincd(cdtrk) .eq. 0) then !Require kincd = 0 * For tracks satisfying a momentum cut, check that the chisquare for * passing through the beam spot passes a cut. Note that we do not * use the returned vertex. if(w(9) .gt. pcut) then !Momentum cut * Test how well track passes through beam vertex call kutl_dcopy(beam_pos, vtx, 3) update = 0 call kvtx_beam(1, ik, update, .TRUE., vtx, Vvtx, chisq, err) chidf = chisq / 2. phideg = ficd(cdtrk) * 180./3.141592654 binphi = min(8, 1 + int(phideg/45.) ) !45 degree phi bins binp = min(10, 1 + int(w(9) / 0.5) ) !0.5 Gev momentum bins * Plot chisq/df vs momentum, phi call hf1(10, real(chidf), 1.) call hf1(10+binp, real(chidf), 1.) call hf1(20+binphi, real(chidf), 1.) * Make list of tracks to use in making vertex if(err.eq.0 .and. chisq.le.chicut) then num_vtx = num_vtx + 1 list_vtx(num_vtx) = ik chisq_in_prim(num_vtx) = chisq chidf_in_prim(num_vtx) = chidf endif endif endif enddo do i=1,3 vtx_prim(i) = bmpos(i,1) enddo call hf1(1, float(num_vtx), 1.) * Plot the errors of the three vertex coordinates do i=1,3 dsigma(i) = sqrt(Vvtx(i,i)) call hf1(101+i, real(dsigma(i)), 1.) enddo * Beam constrained primary vertex. No tracks are updated. if(num_vtx .lt. 2) goto 9999 update = 0 call kvtx_beam(num_vtx, list_vtx, update, .TRUE., vtx_prim, * vvtx_prim, chisq_prim, err) chidf_prim = chisq_prim / (2*num_vtx) num = min(10, num_vtx) call hf1(110, real(chidf_prim), 1.) call hf1(110+num, real(chidf_prim), 1.) * Check chisq of each track relative to the primary vertex just found * by forcing it to pass through the primary vertex. No tracks are * updated. This will get more sophisiticated with time. do it=1,ntrk ik = list(it) cdtrk = kget_track_index(ik) call kget_track_param(ik, w) if(kincd(cdtrk) .eq. 0) then if(w(9) .gt. pcut) then update = 0 call kvtx_fixed(1, ik, update, vtx_prim, chisq, err) chidf = chisq / 2. phideg = ficd(cdtrk) * 180./3.141592654 binphi = min(8, 1 + int(phideg/45.) ) binp = min(10,1 + int(w(9) / 0.5) ) call hf1(210, real(chidf), 1.) call hf1(210+binp, real(chidf), 1.) call hf1(220+binphi, real(chidf), 1.) endif endif enddo if(err) goto 9999 * Normal exit error = 0 return * Not enough tracks 9999 error = 1 return end