************************************************************ subroutine readmipparam include '../common/commontracker.f' include '../common/calib.f' character*60 fname_param 201 format('trk-LADDER',i1,'-mip.dat') do ilad=1,nladders_view write(fname_param,201)ilad print *,'Opening file: ',fname_param open(10, $ FILE='./bin-aux/'//fname_param(1:LNBLNK(fname_param)) $ ,STATUS='UNKNOWN' $ ,IOSTAT=iostat $ ) if(iostat.ne.0)then print*,'READMIPPARAM: *** Error in opening file ***' return endif do iv=1,nviews read(10,* $ ,IOSTAT=iostat $ )pip, $ mip(int(pip),ilad) c print*,ilad,iv,pip,mip(int(pip),ilad) enddo close(10) enddo return end *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** subroutine readchargeparam include '../common/commontracker.f' include '../common/calib.f' character*60 fname_param 201 format('charge-l',i1,'.dat') do ilad=1,nladders_view write(fname_param,201)ilad print *,'Opening file: ',fname_param open(10, $ FILE='./bin-aux/'//fname_param(1:LNBLNK(fname_param)) $ ,STATUS='UNKNOWN' $ ,IOSTAT=iostat $ ) if(iostat.ne.0)then print*,'READCHARGEPARAM: *** Error in opening file ***' return endif do ip=1,nplanes read(10,* $ ,IOSTAT=iostat $ )pip, $ kch(ip,ilad),cch(ip,ilad),sch(ip,ilad) c print*,ilad,ip,pip,kch(ip,ilad), c $ cch(ip,ilad),sch(ip,ilad) enddo close(10) enddo return end *** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** subroutine readetaparam * ----------------------------------------- * read eta2,3,4 calibration parameters * and fill variables: * * eta2(netabin,nladders_view,nviews) * eta3(2*netabin,nladders_view,nviews) * eta4(2*netabin,nladders_view,nviews) * include '../common/commontracker.f' include '../common/calib.f' character*40 fname_binning character*40 fname_param c character*120 cmd1 c character*120 cmd2 ******retrieve ANGULAR BINNING info fname_binning='binning.dat' print *,'Opening file: ',fname_binning open(10, $ FILE='./bin-aux/'//fname_binning(1:LNBLNK(fname_binning)) $ ,STATUS='UNKNOWN' $ ,IOSTAT=iostat $ ) if(iostat.ne.0)then print*,'READETAPARAM: *** Error in opening file ***' return endif print*,'---- ANGULAR BINNING ----' print*,'Bin - angL - angR' 101 format(i2,' ',f6.2,' ',f6.2) do ibin=1,nangmax read(10,* $ ,IOSTAT=iostat $ )xnn,angL(ibin),angR(ibin) if(iostat.ne.0)goto 1000 write(*,101)int(xnn),angL(ibin),angR(ibin) enddo 1000 nangbin=int(xnn) close(10) print*,'-------------------------' do ieta=2,4 !loop on eta 2,3,4 ******retrieve correction parameters 200 format(' Opening eta',i1,' files...') write(*,200)ieta 201 format('eta',i1,'-bin',i1,'-l',i1,'.dat') 202 format('eta',i1,'-bin',i2,'-l',i1,'.dat') do iang=1,nangbin do ilad=1,nladders_view if(iang.lt.10)write(fname_param,201)ieta,iang,ilad if(iang.ge.10)write(fname_param,202)ieta,iang,ilad c print *,'Opening file: ',fname_param open(10, $ FILE='./bin-aux/'//fname_param(1:LNBLNK(fname_param)) $ ,STATUS='UNKNOWN' $ ,IOSTAT=iostat $ ) if(iostat.ne.0)then print*,'READETAPARAM: *** Error in opening file ***' return endif do ival=1,netavalmax if(ieta.eq.2)read(10,* $ ,IOSTAT=iostat $ ) $ eta2(ival,iang), $ (feta2(ival,iv,ilad,iang),iv=1,nviews) if(ieta.eq.3)read(10,* $ ,IOSTAT=iostat $ ) $ eta3(ival,iang), $ (feta3(ival,iv,ilad,iang),iv=1,nviews) if(ieta.eq.4)read(10,* $ ,IOSTAT=iostat $ ) $ eta4(ival,iang), $ (feta4(ival,iv,ilad,iang),iv=1,nviews) if(iostat.ne.0)then netaval=ival-1 c$$$ if(eta2(1,iang).ne.-eta2(netaval,iang)) c$$$ $ print*,'**** ERROR on parameters !!! ****' goto 2000 endif enddo 2000 close(10) * print*,'... done' enddo enddo enddo !end loop on eta 2,3,4 return end ************************************************************ ************************************************************ ************************************************************ ************************************************************ * * This routine provides the coordinates (in cm) in the PAMELA reference system: * - of the point associated with a COUPLE ---> (xPAM,yPAM,zPAM) * - of the extremes of the segment * associated with a SINGLET ---------------> (xPAM_A,yPAM_A,zPAM_A) * ---> (xPAM_B,yPAM_B,zPAM_B) * * It also assigns the spatial resolution to the evaluated coordinates, * as a function (in principle) of the multiplicity, the angle, the PFA etc... * * * To call the routine you must pass the arguments: * icx - ID of cluster x * icy - ID of cluster y * sensor - sensor (1,2) * PFAx - Position Finding Algorithm in x (COG2,ETA2,...) * PFAy - Position Finding Algorithm in y (COG2,ETA2,...) * angx - Projected angle in x * angy - Projected angle in y * * --------- COUPLES ------------------------------------------------------- * The couple defines a point in the space. * The coordinates of the point are evaluated as follows: * 1 - the corrected coordinates relative to the sensor are evaluated * according to the chosen PFA --> (xi,yi,0) * 2 - coordinates are rotated and traslated, according to the aligmnet * parameters, and expressed in the reference system of the mechanical * sensor --> (xrt,yrt,zrt) * 3 - coordinates are finally converted to the PAMELA reference system * --> (xPAM,yPAM,zPAM) * * --------- SINGLETS ------------------------------------------------------- * Since a coordinate is missing, the singlet defines not a point * in the space but a segment AB (parallel to the strips). * In this case the routine returns the coordinates in the PAMELA reference * system of the two extremes A and B of the segment: * --> (xPAM_A,yPAM_A,zPAM_A) * --> (xPAM_B,yPAM_B,zPAM_B) * * ========================================================== * * The output of the routine is stored in the commons: * * double precision xPAM,yPAM,zPAM * common/coord_xyz_PAM/xPAM,yPAM,zPAM * * double precision xPAM_A,yPAM_A,zPAM_A * double precision xPAM_B,yPAM_B,zPAM_B * common/coord_AB_PAM/xPAM_A,yPAM_A,zPAM_A,xPAM_B,yPAM_B,zPAM_B * * double precision resxPAM,resyPAM * common/resolution_PAM/resxPAM,resyPAM * * (in file ../common/common_xyzPAM.f) * * subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,angx,angy) c***************************************************** c 07/10/2005 modified by elena vannuccini --> (1) c 01/02/2006 modified by elena vannuccini --> (2) c 02/02/2006 modified by Elena Vannuccini --> (3) c (implemented new p.f.a.) c 03/02/2006 modified by Elena Vannuccini --> (4) c (implemented variable resolution) c***************************************************** include '../common/commontracker.f' include '../common/calib.f' include '../common/level1.f' include '../common/common_align.f' include '../common/common_mech.f' include '../common/common_xyzPAM.f' include '../common/common_resxy.f' logical DEBUG common/dbg/DEBUG integer icx,icy !X-Y cluster ID integer sensor integer viewx,viewy character*4 PFAx,PFAy !PFA to be used real angx,angy !X-Y angle real stripx,stripy double precision xrt,yrt,zrt double precision xrt_A,yrt_A,zrt_A double precision xrt_B,yrt_B,zrt_B c double precision xi,yi,zi c double precision xi_A,yi_A,zi_A c double precision xi_B,yi_B,zi_B parameter (ndivx=30) resxPAM = 0 resyPAM = 0 xPAM = 0. yPAM = 0. zPAM = 0. xPAM_A = 0. yPAM_A = 0. zPAM_A = 0. xPAM_B = 0. yPAM_B = 0. zPAM_B = 0. * ----------------- * CLUSTER X * ----------------- if(icx.ne.0)then viewx = VIEW(icx) nldx = nld(MAXS(icx),VIEW(icx)) nplx = npl(VIEW(icx)) resxPAM = RESXAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!! stripx = float(MAXS(icx)) if(PFAx.eq.'COG1')then !(1) stripx = stripx !(1) resxPAM = resxPAM !(1) elseif(PFAx.eq.'COG2')then stripx = stripx + cog(2,icx) resxPAM = resxPAM*fbad_cog(2,icx) elseif(PFAx.eq.'ETA2')then c cog2 = cog(2,icx) c etacorr = pfa_eta2(cog2,viewx,nldx,angx) c stripx = stripx + etacorr stripx = stripx + pfa_eta2(icx,angx) !(3) resxPAM = risx_eta2(angx) ! (4) if(DEBUG.and.fbad_cog(2,icx).ne.1) $ print*,'BAD icx >>> ',viewx,fbad_cog(2,icx) resxPAM = resxPAM*fbad_cog(2,icx) elseif(PFAx.eq.'ETA3')then !(3) stripx = stripx + pfa_eta3(icx,angx) !(3) resxPAM = risx_eta3(angx) ! (4) if(DEBUG.and.fbad_cog(3,icx).ne.1) !(3) $ print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)!(3) resxPAM = resxPAM*fbad_cog(3,icx) !(3) elseif(PFAx.eq.'ETA4')then !(3) stripx = stripx + pfa_eta4(icx,angx) !(3) resxPAM = risx_eta4(angx) ! (4) if(DEBUG.and.fbad_cog(4,icx).ne.1) !(3) $ print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)!(3) resxPAM = resxPAM*fbad_cog(4,icx) !(3) elseif(PFAx.eq.'ETA')then !(3) stripx = stripx + pfa_eta(icx,angx) !(3) resxPAM = ris_eta(icx,angx) ! (4) if(DEBUG.and.fbad_cog(2,icx).ne.1) !(3) $ print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)!(3) c resxPAM = resxPAM*fbad_cog(2,icx) !(3)TEMPORANEO resxPAM = resxPAM*fbad_eta(icx,angx) !(3)(4) elseif(PFAx.eq.'COG')then !(2) stripx = stripx + cog(0,icx) !(2) resxPAM = risx_cog(angx) ! (4) resxPAM = resxPAM*fbad_cog(0,icx)!(2) else print*,'*** Non valid p.f.a. (x) --> ',PFAx endif endif * ----------------- * CLUSTER Y * ----------------- if(icy.ne.0)then viewy = VIEW(icy) nldy = nld(MAXS(icy),VIEW(icy)) nply = npl(VIEW(icy)) resyPAM = RESYAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!! if(icx.ne.0.and.(nply.ne.nplx.or.nldy.ne.nldx))then print*,'xyz_PAM ***ERROR*** invalid cluster couple!!! ' $ ,icx,icy goto 100 endif stripy = float(MAXS(icy)) if(PFAy.eq.'COG1')then !(1) stripy = stripy !(1) resyPAM = resyPAM !(1) elseif(PFAy.eq.'COG2')then stripy = stripy + cog(2,icy) resyPAM = resyPAM*fbad_cog(2,icy) elseif(PFAy.eq.'ETA2')then c cog2 = cog(2,icy) c etacorr = pfa_eta2(cog2,viewy,nldy,angy) c stripy = stripy + etacorr stripy = stripy + pfa_eta2(icy,angy) !(3) resyPAM = risy_eta2(angy) ! (4) resyPAM = resyPAM*fbad_cog(2,icy) if(DEBUG.and.fbad_cog(2,icy).ne.1) $ print*,'BAD icy >>> ',viewy,fbad_cog(2,icy) elseif(PFAy.eq.'ETA3')then !(3) stripy = stripy + pfa_eta3(icy,angy) !(3) resyPAM = resyPAM*fbad_cog(3,icy) !(3) if(DEBUG.and.fbad_cog(3,icy).ne.1) !(3) $ print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)!(3) elseif(PFAy.eq.'ETA4')then !(3) stripy = stripy + pfa_eta4(icy,angy) !(3) resyPAM = resyPAM*fbad_cog(4,icy) !(3) if(DEBUG.and.fbad_cog(4,icy).ne.1) !(3) $ print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)!(3) elseif(PFAy.eq.'ETA')then !(3) stripy = stripy + pfa_eta(icy,angy) !(3) resyPAM = ris_eta(icy,angy) ! (4) c resyPAM = resyPAM*fbad_cog(2,icy) !(3)TEMPORANEO resyPAM = resyPAM*fbad_eta(icy,angy) ! (4) if(DEBUG.and.fbad_cog(2,icy).ne.1) !(3) $ print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)!(3) elseif(PFAy.eq.'COG')then stripy = stripy + cog(0,icy) resyPAM = risy_cog(angy) ! (4) c resyPAM = ris_eta(icy,angy) ! (4) resyPAM = resyPAM*fbad_cog(0,icy) else print*,'*** Non valid p.f.a. (x) --> ',PFAx endif endif c=========================================================== C COUPLE C=========================================================== if(icx.ne.0.and.icy.ne.0)then c------------------------------------------------------------------------ c (xi,yi,zi) = mechanical coordinates in the silicon sensor frame c------------------------------------------------------------------------ xi = acoordsi(stripx,viewx) yi = acoordsi(stripy,viewy) zi = 0. c------------------------------------------------------------------------ c (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame c------------------------------------------------------------------------ c N.B. I convert angles from microradiants to radiants xrt = xi $ - omega(nplx,nldx,sensor)*yi $ + gamma(nplx,nldx,sensor)*zi $ + dx(nplx,nldx,sensor) yrt = omega(nplx,nldx,sensor)*xi $ + yi $ - beta(nplx,nldx,sensor)*zi $ + dy(nplx,nldx,sensor) zrt = -gamma(nplx,nldx,sensor)*xi $ + beta(nplx,nldx,sensor)*yi $ + zi $ + dz(nplx,nldx,sensor) c xrt = xi c yrt = yi c zrt = zi c------------------------------------------------------------------------ c (xPAM,yPAM,zPAM) = measured coordinates (in cm) c in PAMELA reference system c------------------------------------------------------------------------ xPAM = dcoord(xrt,viewx,nldx,sensor) / 1.d4 yPAM = dcoord(yrt,viewy,nldy,sensor) / 1.d4 zPAM = ( zrt + z_mech_sensor(nplx,nldx,sensor)*1000. ) / 1.d4 xPAM_A = 0. yPAM_A = 0. zPAM_A = 0. xPAM_B = 0. yPAM_B = 0. zPAM_B = 0. elseif( $ (icx.ne.0.and.icy.eq.0).or. $ (icx.eq.0.and.icy.ne.0).or. $ .false. $ )then c------------------------------------------------------------------------ c (xi,yi,zi) = mechanical coordinates in the silicon sensor frame c------------------------------------------------------------------------ if(icy.ne.0)then c=========================================================== C Y-SINGLET C=========================================================== nplx = nply nldx = nldy viewx = viewy + 1 yi = acoordsi(stripy,viewy) xi_A = edgeY_d - SiDimX/2 yi_A = yi zi_A = 0. xi_B = SiDimX/2 - edgeY_u yi_B = yi zi_B = 0. c print*,'Y-cl ',icy,stripy,' --> ',yi c print*,xi_A,' <--> ',xi_B elseif(icx.ne.0)then c=========================================================== C X-SINGLET C=========================================================== nply = nplx nldy = nldx viewy = viewx - 1 xi = acoordsi(stripx,viewx) xi_A = xi yi_A = edgeX_d - SiDimY/2 zi_A = 0. xi_B = xi yi_B = SiDimY/2 - edgeX_u zi_B = 0. if(viewy.eq.11)then yi = yi_A yi_A = yi_B yi_B = yi endif c print*,'X-cl ',icx,stripx,' --> ',xi c print*,yi_A,' <--> ',yi_B else print *,'routine xyz_PAM ---> not properly used !!!' print *,'icx = ',icx print *,'icy = ',icy goto 100 endif c------------------------------------------------------------------------ c (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame c------------------------------------------------------------------------ c N.B. I convert angles from microradiants to radiants xrt_A = xi_A $ - omega(nplx,nldx,sensor)*yi_A $ + gamma(nplx,nldx,sensor)*zi_A $ + dx(nplx,nldx,sensor) yrt_A = omega(nplx,nldx,sensor)*xi_A $ + yi_A $ - beta(nplx,nldx,sensor)*zi_A $ + dy(nplx,nldx,sensor) zrt_A = -gamma(nplx,nldx,sensor)*xi_A $ + beta(nplx,nldx,sensor)*yi_A $ + zi_A $ + dz(nplx,nldx,sensor) xrt_B = xi_B $ - omega(nplx,nldx,sensor)*yi_B $ + gamma(nplx,nldx,sensor)*zi_B $ + dx(nplx,nldx,sensor) yrt_B = omega(nplx,nldx,sensor)*xi_B $ + yi_B $ - beta(nplx,nldx,sensor)*zi_B $ + dy(nplx,nldx,sensor) zrt_B = -gamma(nplx,nldx,sensor)*xi_B $ + beta(nplx,nldx,sensor)*yi_B $ + zi_B $ + dz(nplx,nldx,sensor) c xrt = xi c yrt = yi c zrt = zi c------------------------------------------------------------------------ c (xPAM,yPAM,zPAM) = measured coordinates (in cm) c in PAMELA reference system c------------------------------------------------------------------------ xPAM = 0. yPAM = 0. zPAM = 0. xPAM_A = dcoord(xrt_A,viewx,nldx,sensor) / 1.d4 yPAM_A = dcoord(yrt_A,viewy,nldy,sensor) / 1.d4 zPAM_A = ( zrt_A + z_mech_sensor(nplx,nldx,sensor)*1000.)/ 1.d4 xPAM_B = dcoord(xrt_B,viewx,nldx,sensor) / 1.d4 yPAM_B = dcoord(yrt_B,viewy,nldy,sensor) / 1.d4 zPAM_B = ( zrt_B + z_mech_sensor(nplx,nldx,sensor)*1000.)/ 1.d4 c print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')' else print *,'routine xyz_PAM ---> not properly used !!!' print *,'icx = ',icx print *,'icy = ',icy endif 100 continue end ******************************************************************************** ******************************************************************************** ******************************************************************************** * * The function distance_to(XP,YP) should be used after * a call to the xyz_PAM routine and it evaluate the * NORMALIZED distance (PROJECTED on the XY plane) between * the point (XP,YP), argument of the function, * and: * * - the point (xPAM,yPAM,zPAM), in the case of a COUPLE * or * - the segment (xPAM_A,yPAM_A,zPAM_A)-(xPAM_B,yPAM_B,zPAM_B), * in the case of a SINGLET. * * ( The routine xyz_PAM fills the common defined in "common_xyzPAM.f", * which stores the coordinates of the couple/singlet ) * ******************************************************************************** real function distance_to(XPP,YPP) include '../common/common_xyzPAM.f' * ----------------------------------- * it computes the normalized distance * ( i.e. distance/resolution ) * ----------------------------------- double precision distance,RE double precision BETA,ALFA,xmi,ymi * ---------------------- if ( + xPAM.eq.0.and. + yPAM.eq.0.and. + zPAM.eq.0.and. + xPAM_A.ne.0.and. + yPAM_A.ne.0.and. + zPAM_A.ne.0.and. + xPAM_B.ne.0.and. + yPAM_B.ne.0.and. + zPAM_B.ne.0.and. + .true.)then * ----------------------- * DISTANCE TO --- SINGLET * ----------------------- if(abs(sngl(xPAM_B-xPAM_A)).lt.abs(sngl(yPAM_B-yPAM_A)))then * |||---------- X CLUSTER BETA = (xPAM_B-xPAM_A)/(yPAM_B-yPAM_A) ALFA = xPAM_A - BETA * yPAM_A ymi = ( YPP + BETA*XPP - BETA*ALFA )/(1+BETA**2) if(ymi.lt.dmin1(yPAM_A,yPAM_B))ymi=dmin1(yPAM_A,yPAM_B) if(ymi.gt.dmax1(yPAM_A,yPAM_B))ymi=dmax1(yPAM_A,yPAM_B) xmi = ALFA + BETA * ymi RE = resxPAM else * |||---------- Y CLUSTER BETA = (yPAM_B-yPAM_A)/(xPAM_B-xPAM_A) ALFA = yPAM_A - BETA * xPAM_A xmi = ( XPP + BETA*YPP - BETA*ALFA )/(1+BETA**2) if(xmi.lt.dmin1(xPAM_A,xPAM_B))xmi=dmin1(xPAM_A,xPAM_B) if(xmi.gt.dmax1(xPAM_A,xPAM_B))xmi=dmax1(xPAM_A,xPAM_B) ymi = ALFA + BETA * xmi RE = resyPAM endif distance= $ ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2 distance=dsqrt(distance) c$$$ print*,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b c$$$ $ ,' --- distance_to --- ',xpp,ypp c$$$ print*,' resolution ',re * ---------------------- elseif( + xPAM.ne.0.and. + yPAM.ne.0.and. + zPAM.ne.0.and. + xPAM_A.eq.0.and. + yPAM_A.eq.0.and. + zPAM_A.eq.0.and. + xPAM_B.eq.0.and. + yPAM_B.eq.0.and. + zPAM_B.eq.0.and. + .true.)then * ---------------------- * DISTANCE TO --- COUPLE * ---------------------- distance= $ ((xPAM-XPP)/resxPAM)**2 $ + $ ((yPAM-YPP)/resyPAM)**2 distance=dsqrt(distance) c$$$ print*,xPAM,yPAM,zPAM c$$$ $ ,' --- distance_to --- ',xpp,ypp c$$$ print*,' resolution ',resxPAM,resyPAM else print* $ ,' function distance_to ---> wrong usage!!!' print*,' xPAM,yPAM,zPAM ',xPAM,yPAM,zPAM print*,' xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b ' $ ,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b endif distance_to = sngl(distance) return end ******************************************************************************** ******************************************************************************** ******************************************************************************** ******************************************************************************** subroutine whichsensor(nplPAM,xPAM,yPAM,ladder,sensor) * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Given the plane (1-6 from BOTTOM to TOP!!) and the (xPAM,yPAM) * coordinates (in the PAMELA reference system), it returns * the ladder and the sensor which the point belongs to. * * The method to assign a point to a sensor consists in * - calculating the sum of the distances between the point * and the sensor edges * - requiring that it is less-equal than (SiDimX+SiDimY) * * NB -- SiDimX and SiDimY are not the dimentions of the SENSITIVE volume * but of the whole silicon sensor * * CONVENTION: * - sensor 1 is the one closest to the hybrid * - ladder 1 is the first to be read out (strips from 1 to 1024) * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * include '../common/commontracker.f' include '../common/common_align.f' integer ladder,sensor,viewx,viewy real c1(4),c2(4),c3(4) data c1/1.,0.,0.,1./ data c2/1.,-1.,-1.,1./ data c3/1.,1.,0.,0./ real*8 yvvv,xvvv double precision xi,yi,zi double precision xrt,yrt,zrt real AA,BB real yvv(4),xvv(4) * tollerance to consider the track inside the sensitive area real ptoll data ptoll/150./ !um external nviewx,nviewy,acoordsi,dcoord nplpt = nplPAM !plane viewx = nviewx(nplpt) viewy = nviewy(nplpt) do il=1,nladders_view do is=1,2 do iv=1,4 !loop on sensor vertexes stripx = (il-c1(iv))*1024 + c1(iv) + c2(iv)*3 stripy = (il-c3(iv))*1024 + c3(iv) c------------------------------------------------------------------------ c (xi,yi,zi) = mechanical coordinates in the silicon sensor frame c------------------------------------------------------------------------ xi = acoordsi(stripx,viewx) yi = acoordsi(stripy,viewy) zi = 0. c------------------------------------------------------------------------ c (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame c------------------------------------------------------------------------ c N.B. I convert angles from microradiants to radiants xrt = xi $ - omega(nplpt,il,is)*yi $ + gamma(nplpt,il,is)*zi $ + dx(nplpt,il,is) yrt = omega(nplpt,il,is)*xi $ + yi $ - beta(nplpt,il,is)*zi $ + dy(nplpt,il,is) zrt = -gamma(nplpt,il,is)*xi $ + beta(nplpt,il,is)*yi $ + zi $ + dz(nplpt,il,is) c------------------------------------------------------------------------ c measured coordinates (in cm) in PAMELA reference system c------------------------------------------------------------------------ yvvv = dcoord(yrt,viewy,il,is) / 1.d4 xvvv = dcoord(xrt,viewx,il,is) / 1.d4 yvv(iv)=sngl(yvvv) xvv(iv)=sngl(xvvv) c print*,'LADDER ',il,' SENSOR ',is,' vertexes >> ' c $ ,iv,xvv(iv),yvv(iv) enddo !end loop on sensor vertexes dtot=0. do iside=1,4,2 !loop on sensor edges X iv1=iside iv2=mod(iside,4)+1 * straight line passing trhough two consecutive vertexes AA = (yvv(iv1)-yvv(iv2))/(xvv(iv1)-xvv(iv2)) BB = yvv(iv1) - AA*xvv(iv1) * point along the straight line closer to the track xoo = (xPAM+AA*yPAM-AA*BB)/(1+AA**2) yoo = AA*xoo + BB * sum of the distances dtot = dtot + $ sqrt((xPAM-xoo)**2+(yPAM-yoo)**2) enddo !end loop on sensor edges do iside=2,4,2 !loop on sensor edges Y iv1=iside iv2=mod(iside,4)+1 * straight line passing trhough two consecutive vertexes AA = (xvv(iv1)-xvv(iv2))/(yvv(iv1)-yvv(iv2)) BB = xvv(iv1) - AA*yvv(iv1) * point along the straight line closer to the track yoo = (yPAM+AA*xPAM-AA*BB)/(1+AA**2) xoo = AA*yoo + BB * sum of the distances dtot = dtot + $ sqrt((xPAM-xoo)**2+(yPAM-yoo)**2) enddo !end loop on sensor edges * half-perimeter of sensitive area Perim = $ SiDimX - edgeX_l - edgeX_r $ +SiDimY - edgeY_l - edgeY_r Perim = (Perim + ptoll)/1.e4 if(dtot.le.Perim)goto 100 enddo enddo ladder = 0 sensor = 0 goto 200 100 continue ladder = il sensor = is 200 return end ************************************************************************* subroutine reverse(v,n,temp) !invert the order of the components of v(n) vector implicit double precision (A-H,O-Z) dimension v(*) dimension temp(*) integer i,n do i=1,n temp(i)=v(n+1-i) enddo do i=1,n v(i)=temp(i) enddo return end ************************************************************************* ************************************************************************* ************************************************************************* ************************************************************************* ************************************************************************* ************************************************************************* ************************************************************************* ************************************************************************* ************************************************************************* ************************************************************************* ************************************************************************* ************************************************************************* ************************************************************************* ************************************************************************* ************************************************************************* integer function ip_cp(id) * * given the couple id, * it returns the plane number * include '../common/commontracker.f' c include '../common/common_analysis.f' include '../common/common_momanhough.f' ip_cp=0 ncpp=0 do ip=1,nplanes ncpp=ncpp+ncp_plane(ip) if(ncpp.ge.abs(id))then ip_cp=ip goto 100 endif enddo 100 continue return end integer function is_cp(id) * * given the couple id, * it returns the sensor number * is_cp=0 if(id.lt.0)is_cp=1 if(id.gt.0)is_cp=2 if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!' return end integer function icp_cp(id) * * given the couple id, * it returns the id number ON THE PLANE * include '../common/commontracker.f' c include '../common/common_analysis.f' include '../common/common_momanhough.f' icp_cp=0 ncpp=0 do ip=1,nplanes ncppold=ncpp ncpp=ncpp+ncp_plane(ip) if(ncpp.ge.abs(id))then icp_cp=abs(id)-ncppold goto 100 endif enddo 100 continue return end integer function id_cp(ip,icp,is) * * given a plane, a couple and the sensor * it returns the absolute couple id * negative if sensor =1 * positive if sensor =2 * include '../common/commontracker.f' c include '../common/calib.f' c include '../common/level1.f' c include '../common/common_analysis.f' include '../common/common_momanhough.f' id_cp=0 if(ip.gt.1)then do i=1,ip-1 id_cp = id_cp + ncp_plane(i) enddo endif id_cp = id_cp + icp if(is.eq.1) id_cp = -id_cp return end ************************************************************************* ************************************************************************* ************************************************************************* ************************************************************************* ************************************************************************* ************************************************************************* subroutine book_debug include '../common/commontracker.f' include '../common/common_momanhough.f' include '../common/common_level2debug.f' character*35 block1,block2,block3!,block4 $ ,block5!,block6 * * * * * * * * * * * * * * * * * * * * * * * * * HOUGH TRANSFORM PARAMETERS call HBOOK2(1003 $ ,'y vs tg thyz' $ ,300,-1.,1. !x $ ,3000,-70.,70.,0.) !y call HBOOK1(1004 $ ,'Dy' $ ,100,0.,2.,0.) call HBOOK1(1005 $ ,'D thyz' $ ,100,0.,.05,0.) * DEBUG ntuple: call HBNT(ntp_level2+1,'LEVEL2',' ') call HBNAME(ntp_level2+1,'EVENT',good2_nt, $ 'GOOD2:L,NEV2:I') 411 format('NDBLT:I::[0,',I5,']') write(block1,411) ndblt_max_nt call HBNAME(ntp_level2+1,'HOUGH YZ',ndblt_nt, $ block1//' $ ,ALFAYZ1(NDBLT):R $ ,ALFAYZ2(NDBLT):R $ ,DB_CLOUD(NDBLT):I $ ') 412 format('NTRPT:I::[0,',I5,']') write(block2,412) ntrpt_max_nt call HBNAME(ntp_level2+1,'HOUGH XZ',NTRPT_nt, $ block2//' $ ,ALFAXZ1(NTRPT):R $ ,ALFAXZ2(NTRPT):R $ ,ALFAXZ3(NTRPT):R $ ,TR_CLOUD(NTRPT):I $ ') 413 format('NCLOUDS_YZ:I::[0,',I4,']') c$$$ 414 format('DB_CLOUD(',I4,'):I') write(block3,413) ncloyz_max c$$$ write(block4,414) ndblt_max_nt call HBNAME(ntp_level2+1,'CLOUD YZ',NCLOUDS_YZ, $ block3//' $ ,ALFAYZ1_AV(NCLOUDS_YZ):R $ ,ALFAYZ2_AV(NCLOUDS_YZ):R $ ,PTCLOUD_YZ(NCLOUDS_YZ):I' c$$$ $ ,'//block4 $ ) 415 format('NCLOUDS_XZ:I::[0,',I4,']') c$$$ 416 format('TR_CLOUD(',I5,'):I') write(block5,415) ncloxz_max c$$$ write(block6,416) ntrpt_max_nt call HBNAME(ntp_level2+1,'CLOUD XZ',NCLOUDS_XZ, $ block5//' $ ,ALFAXZ1_AV(NCLOUDS_XZ):R $ ,ALFAXZ2_AV(NCLOUDS_XZ):R $ ,ALFAXZ3_AV(NCLOUDS_XZ):R $ ,PTCLOUD_XZ(NCLOUDS_XZ):I' c$$$ $ ,'//block6 $ ) return end ***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...*** * * * * * * * * * ***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...*** subroutine book_level2 c***************************************************** cccccc 11/9/2005 modified by david fedele cccccc 07/10/2005 modified by elena vannuccini --> (2) c***************************************************** include '../common/commontracker.f' include '../common/common_momanhough.f' include '../common/level2.f' character*35 block1,block2 c print*,'__________ booking LEVEL2 n-tuple __________' * LEVEL1 ntuple: call HBNT(ntp_level2,'LEVEL2',' ') c***************************************************** cccccc 11/9/2005 modified by david fedele c call HBNAME(ntp_level2,'EVENT',good2,'GOOD2:L,NEV2:I') cccccc 06/10/2005 modified by elena vannuccini c call HBNAME(ntp_level2,'GENERAL',good2,'GOOD2:L,NEV2:I c $ ,WHIC_CALIB:I,SWCODE:I') call HBNAME(ntp_level2,'GENERAL',good2,'GOOD2:L,NEV2:I $ ,WHICH_CALIB:I,SWCODE:I,CRC(12):L') c********************************************************* # ifndef TEST2003 call HBNAME(ntp_level2,'CPU',pkt_type $ ,'PKT_TYPE:I::[0,50] $ ,PKT_NUM:I $ ,OBT:I'// c******************************************************** cccccc 11/9/2005 modified by david fedele c $ ,WHICH_CALIB:I::[0,50]') $ ',CPU_CRC:L') c******************************************************** # endif 417 format('NTRK:I::[0,',I4,']') 418 format(',IMAGE(NTRK):I::[0,',I4,']') write(block1,417)NTRKMAX write(block2,418)NTRKMAX call HBNAME(ntp_level2,'TRACKS',NTRK, $ block1// $ block2//' $ ,XM(6,NTRK):R $ ,YM(6,NTRK):R $ ,ZM(6,NTRK):R $ ,RESX(6,NTRK):R $ ,RESY(6,NTRK):R $ ,AL(5,NTRK):R $ ,COVAL(5,5,NTRK):R $ ,CHI2(NTRK):R $ ,XGOOD(6,NTRK):I::[0,1] $ ,YGOOD(6,NTRK):I::[0,1] $ ,XV(6,NTRK):R $ ,YV(6,NTRK):R $ ,ZV(6,NTRK):R $ ,AXV(6,NTRK):R $ ,AYV(6,NTRK):R'// c***************************************************** cccccc 11/9/2005 modified by david fedele c $ ,DEDXP(6,NTRK):R'// c $ ') $ ',DEDX_X(6,NTRK):R $ ,DEDX_Y(6,NTRK):R'// c**************************************************** cccccc 06/10/2005 modified by elena vannuccini c $ ,CRC(12):L $ ',BdL(NTRK):R' $ ) c**************************************************** call HBNAME(ntp_level2,'SINGLETX',nclsx, c***************************************************** cccccc 11/9/2005 modified by david fedele c $ 'NCLSX(6):I,NCLSY(6):I') $ 'NCLSX:I::[0,500],PLANEX(NCLSX):I $ ,XS(2,NCLSX):R,SGNLXS(NCLSX):R') !(2) c $ ,XS(NCLSX):R,SGNLXS(NCLSX):R') !(2) call HBNAME(ntp_level2,'SINGLETY',nclsy, $ 'NCLSY:I::[0,500],PLANEY(NCLSY):I $ ,YS(2,NCLSY):R,SGNLYS(NCLSY):R') !(2) c $ ,YS(NCLSY):R,SGNLYS(NCLSY):R') !(2) return end * **************************************************** subroutine init_level2 c***************************************************** c 07/10/2005 modified by elena vannuccini --> (1) c***************************************************** include '../common/commontracker.f' include '../common/common_momanhough.f' include '../common/level2.f' include '../common/level1.f' good2 = .false. nev2 = nev1 # ifndef TEST2003 c***************************************************** cccccc 11/9/2005 modified by david fedele c pkt_type = pkt_type1 c pkt_num = pkt_num1 c obt = obt1 c which_calib = which_calib1 swcode = 302 which_calib = which_calib1 pkt_type = pkt_type1 pkt_num = pkt_num1 obt = obt1 cpu_crc = cpu_crc1 do iv=1,12 crc(iv)=crc1(iv) enddo # endif c***************************************************** NTRK = 0 do it=1,NTRKMAX!NTRACKSMAX IMAGE(IT)=0 CHI2_NT(IT) = -100000. BdL(IT) = 0. do ip=1,nplanes XM_nt(IP,IT) = 0 YM_nt(IP,IT) = 0 ZM_nt(IP,IT) = 0 RESX_nt(IP,IT) = 0 RESY_nt(IP,IT) = 0 XGOOD_nt(IP,IT) = 0 YGOOD_nt(IP,IT) = 0 c***************************************************** cccccc 11/9/2005 modified by david fedele DEDX_X(IP,IT) = 0 DEDX_Y(IP,IT) = 0 c****************************************************** enddo do ipa=1,5 AL_nt(IPA,IT) = 0 do ipaa=1,5 coval(ipa,ipaa,IT)=0 enddo enddo enddo c***************************************************** cccccc 11/9/2005 modified by david fedele nclsx=0 nclsy=0 do ip=1,NSINGMAX planex(ip)=0 c xs(ip)=0 xs(1,ip)=0 xs(2,ip)=0 sgnlxs(ip)=0 planey(ip)=0 c ys(ip)=0 ys(1,ip)=0 ys(2,ip)=0 sgnlys(ip)=0 enddo c******************************************************* end subroutine init_level2_debug c***************************************************** c 01/12/2005 createded by elena vannuccini c***************************************************** include '../common/commontracker.f' include '../common/common_momanhough.f' include '../common/common_level2debug.f' include '../common/level2.f' c include '../common/level1.f' good2_nt = .false. nev2_nt = nev2 ntrpt_nt=0 ndblt_nt=0 NCLOUDS_XZ_nt=0 NCLOUDS_YZ_nt=0 do idb=1,ndblt_max_nt db_cloud_nt(idb)=0 enddo do itr=1,ntrpl_max_nt tr_cloud_nt(itr)=0 enddo end ************************************************************ * * * * * * * ************************************************************ subroutine fill_level2_tracks(ntr) * ------------------------------------------------------- * This routine fills the ntr-th element of the variables * inside the level2_tracks common, which correspond * to the ntr-th track info. * ------------------------------------------------------- c***************************************************** cccccc 27/09/2005 modified by david fedele ---> (1) c to normalize al(3) and al(4) cccccc 07/10/2005 modified by elena vannuccini ---> (2) c to fill the new variables cccccc 12/10/2005 modified by elena vannuccini ---> (3) c to fill the BdL variable cccccc 20/10/2005 modified by elena vannuccini ---> (4) c bug in angular track-parameters c (tracking is wrong) c***************************************************** include '../common/commontracker.f' include '../common/level2.f' include '../common/common_mini_2.f' c real gamma,delta !(1) real sinth,phi,pig !(4) pig=acos(-1.) good2=.true. chi2_nt(ntr) = sngl(chi2) c print*,chi2_nt(ntr) c***************************************************** cccccc (4) c***************************************************** c$$$ delta=al(4) !(1) c$$$ gamma=al(3) !(1) c$$$ if (cos(delta).gt.0) then !(1) c$$$ do idiv=1,100 !(1) c$$$ if (delta.gt.(PIGR/2)) then !(1) c$$$ delta=mod(al(4),(idiv*2*PIGR)) !(1) c$$$ endif !(1) c$$$ if((delta.lt.(PIGR/2)).and.(delta.gt.-(PIGR/2))) then !(1) c$$$ al(4)=delta !(1) c$$$ goto 42 !(1) c$$$ endif !(1) c$$$ enddo !(1) c$$$ elseif(cos(delta).lt.0) then !(1) c$$$ al(3)=-gamma !(1) c$$$ do idiv=1,100 !(1) c$$$ if (delta.gt.(PIGR/2)) then !(1) c$$$ delta=mod(al(4),(idiv*3*PIGR)) !(1) c$$$ endif !(1) c$$$ if((delta.lt.(3*PIGR/2)).and.(delta.gt.(PIGR/2)))then !(1) c$$$ al(4)=delta !(1) c$$$ goto 42 !(1) c$$$ endif !(1) c$$$ enddo !(1) c$$$ endif !(1) c$$$ c$$$ 42 continue !(1) ***************************************************** phi = al(4) !(4) sinth = al(3) !(4) if(sinth.lt.0)then !(4) sinth = -sinth !(4) phi = phi + pig !(4) endif !(4) npig = aint(phi/(2*pig)) !(4) phi = phi - npig*2*pig !(4) if(phi.lt.0) !(4) $ phi = phi + 2*pig !(4) al(4) = phi !(4) al(3) = sinth !(4) ***************************************************** do i=1,5 al_nt(i,ntr) = sngl(al(i)) do j=1,5 coval(i,j,ntr) = sngl(cov(i,j)) enddo c print*,al_nt(i,ntr) enddo do ip=1,nplanes ! loop on planes xgood_nt(ip,ntr) = int(xgood(ip)) ygood_nt(ip,ntr) = int(ygood(ip)) xm_nt(ip,ntr) = sngl(xm(ip)) ym_nt(ip,ntr) = sngl(ym(ip)) zm_nt(ip,ntr) = sngl(zm(ip)) RESX_nt(IP,ntr) = sngl(resx(ip)) RESY_nt(IP,ntr) = sngl(resy(ip)) xv_nt(ip,ntr) = sngl(xv(ip)) yv_nt(ip,ntr) = sngl(yv(ip)) zv_nt(ip,ntr) = sngl(zv(ip)) axv_nt(ip,ntr) = sngl(axv(ip)) ayv_nt(ip,ntr) = sngl(ayv(ip)) c dedxp(ip,ntr) = sngl(dedxtrk(ip)) !(1) dedx_x(ip,ntr) = sngl(dedxtrk_x(ip)) !(2) dedx_y(ip,ntr) = sngl(dedxtrk_y(ip)) !(2) enddo call CalcBdL(100,xxxx,IFAIL) if(ifps(xxxx).eq.1)BdL(ntr) = xxxx end subroutine fill_level2_siglets c***************************************************** c 07/10/2005 created by elena vannuccini c 31/01/2006 modified by elena vannuccini * to convert adc to mip --> (2) c***************************************************** * ------------------------------------------------------- * This routine fills the elements of the variables * inside the level2_singletsx and level2_singletsy commons, * which store info on clusters outside the tracks * ------------------------------------------------------- include '../common/commontracker.f' include '../common/level1.f' include '../common/level2.f' include '../common/calib.f' include '../common/common_momanhough.f' include '../common/common_xyzPAM.f' * count #cluster per plane not associated to any track good2=.true. nclsx = 0 nclsy = 0 do icl=1,nclstr1 if(cl_used(icl).eq.0)then !cluster not included in any track ip=nplanes-npl(VIEW(icl))+1 if(mod(VIEW(icl),2).eq.0)then !=== X views nclsx = nclsx + 1 planex(nclsx) = ip sgnlxs(nclsx) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2) do is=1,2 c call xyz_PAM(icl,0,is,'COG1',' ',0.,0.) call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.) xs(is,nclsx) = (xPAM_A+xPAM_B)/2 enddo else !=== Y views nclsy = nclsy + 1 planey(nclsy) = ip sgnlys(nclsy) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2) do is=1,2 c call xyz_PAM(0,icl,is,' ','COG1',0.,0.) call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.) ys(is,nclsy) = (yPAM_A+yPAM_B)/2 enddo endif endif c print*,icl,cl_used(icl),cl_good(icl),ip,VIEW(icl)!nclsx(ip),nclsy(ip) enddo end subroutine fill_level2_clouds c***************************************************** c 29/11/2005 created by elena vannuccini c***************************************************** * ------------------------------------------------------- * This routine fills the variables related to the hough * transform, for the debig n-tuple * ------------------------------------------------------- include '../common/commontracker.f' include '../common/common_momanhough.f' include '../common/common_level2debug.f' include '../common/level2.f' good2_nt=.true.!good2 c nev2_nt=nev2 if(.false. $ .or.ntrpt.gt.ntrpt_max_nt $ .or.ndblt.gt.ndblt_max_nt $ .or.NCLOUDS_XZ.gt.ncloxz_max $ .or.NCLOUDS_yZ.gt.ncloyz_max $ )then good2_nt=.false. ntrpt_nt=0 ndblt_nt=0 NCLOUDS_XZ_nt=0 NCLOUDS_YZ_nt=0 else ndblt_nt=ndblt ntrpt_nt=ntrpt if(ndblt.ne.0)then do id=1,ndblt alfayz1_nt(id)=alfayz1(id) !Y0 alfayz2_nt(id)=alfayz2(id) !tg theta-yz c db_cloud_nt(id)=db_cloud(id) enddo endif if(ndblt.ne.0)then do it=1,ntrpt alfaxz1_nt(it)=alfaxz1(it) !X0 alfaxz2_nt(it)=alfaxz2(it) !tg theta-xz alfaxz3_nt(it)=alfaxz3(it) !1/r c tr_cloud_nt(it)=tr_cloud(it) enddo endif nclouds_yz_nt=nclouds_yz nclouds_xz_nt=nclouds_xz if(nclouds_yz.ne.0)then nnn=0 do iyz=1,nclouds_yz ptcloud_yz_nt(iyz)=ptcloud_yz(iyz) alfayz1_av_nt(iyz)=alfayz1_av(iyz) alfayz2_av_nt(iyz)=alfayz2_av(iyz) nnn=nnn+ptcloud_yz(iyz) enddo do ipt=1,nnn db_cloud_nt(ipt)=db_cloud(ipt) enddo c print*,'#### ntupla #### ptcloud_yz ' c $ ,(ptcloud_yz(i),i=1,nclouds_yz) c print*,'#### ntupla #### db_cloud ' c $ ,(db_cloud(i),i=1,nnn) endif if(nclouds_xz.ne.0)then nnn=0 do ixz=1,nclouds_xz ptcloud_xz_nt(ixz)=ptcloud_xz(ixz) alfaxz1_av_nt(ixz)=alfaxz1_av(ixz) alfaxz2_av_nt(ixz)=alfaxz2_av(ixz) alfaxz3_av_nt(ixz)=alfaxz3_av(ixz) nnn=nnn+ptcloud_xz(ixz) enddo do ipt=1,nnn tr_cloud_nt(ipt)=tr_cloud(ipt) enddo c print*,'#### ntupla #### ptcloud_xz ' c $ ,(ptcloud_xz(i),i=1,nclouds_xz) c print*,'#### ntupla #### tr_cloud ' c $ ,(tr_cloud(i),i=1,nnn) endif endif end *************************************************** * * * * * * * * * * * * ************************************************** subroutine cl_to_couples(iflag) include '../common/commontracker.f' include '../common/common_momanhough.f' include '../common/momanhough_init.f' include '../common/calib.f' include '../common/level1.f' logical DEBUG common/dbg/DEBUG * output flag * -------------- * 0 = good event * 1 = bad event * -------------- integer iflag integer badseed,badcl * init variables ncp_tot=0 do ip=1,nplanes do ico=1,ncouplemax clx(ip,ico)=0 cly(ip,ico)=0 enddo ncp_plane(ip)=0 do icl=1,nclstrmax_level2 cls(ip,icl)=1 enddo ncls(ip)=0 enddo do icl=1,nclstrmax_level2 cl_single(icl)=1 cl_good(icl)=0 enddo * start association ncouples=0 do icx=1,nclstr1 !loop on cluster (X) if(mod(VIEW(icx),2).eq.1)goto 10 * ---------------------------------------------------- * cut on charge (X VIEW) if(dedx(icx).lt.dedx_x_min)then cl_single(icx)=0 goto 10 endif * cut BAD (X VIEW) badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx))) ifirst=INDSTART(icx) if(icx.ne.nclstr1) then ilast=INDSTART(icx+1)-1 else ilast=TOTCLLENGTH endif badcl=badseed do igood=-ngoodstr,ngoodstr ibad=1 if((INDMAX(icx)+igood).gt.ifirst.and. $ (INDMAX(icx)+igood).lt.ilast.and. $ .true.)then ibad=BAD(VIEW(icx), $ nvk(MAXS(icx)+igood), $ nst(MAXS(icx)+igood)) endif badcl=badcl*ibad enddo c if(badcl.eq.0)then c cl_single(icx)=0 c goto 10 c endif * ---------------------------------------------------- cl_good(icx)=1 nplx=npl(VIEW(icx)) nldx=nld(MAXS(icx),VIEW(icx)) do icy=1,nclstr1 !loop on cluster (Y) if(mod(VIEW(icy),2).eq.0)goto 20 * ---------------------------------------------------- * cut on charge (Y VIEW) if(dedx(icy).lt.dedx_y_min)then cl_single(icy)=0 goto 20 endif * cut BAD (Y VIEW) badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy))) ifirst=INDSTART(icy) if(icy.ne.nclstr1) then ilast=INDSTART(icy+1)-1 else ilast=TOTCLLENGTH endif badcl=badseed do igood=-ngoodstr,ngoodstr ibad=1 if((INDMAX(icy)+igood).gt.ifirst.and. $ (INDMAX(icy)+igood).lt.ilast.and. $ .true.) $ ibad=BAD(VIEW(icy), $ nvk(MAXS(icy)+igood), $ nst(MAXS(icy)+igood)) badcl=badcl*ibad enddo c if(badcl.eq.0)then c cl_single(icy)=0 c goto 20 c endif * ---------------------------------------------------- cl_good(icy)=1 nply=npl(VIEW(icy)) nldy=nld(MAXS(icy),VIEW(icy)) * ---------------------------------------------- * CONDITION TO FORM A COUPLE * ---------------------------------------------- * geometrical consistency (same plane and ladder) if(nply.eq.nplx.and.nldy.eq.nldx)then * charge correlation ddd=(dedx(icy) $ -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx)) ddd=ddd/sqrt(kch(nplx,nldx)**2+1) cut=chcut*sch(nplx,nldx) if(abs(ddd).gt.cut)goto 20 !charge not consistent * ------------------> COUPLE <------------------ * check to do not overflow vector dimentions if(ncp_plane(nplx).gt.ncouplemax)then if(DEBUG)print*, $ ' ** warning ** number of identified'// $ ' couples on plane ',nplx, $ ' exceeds vector dimention'// $ ' ( ',ncouplemax,' )' c good2=.false. c goto 880 !fill ntp and go to next event iflag=1 return endif if(ncp_plane(nplx).eq.ncouplemax)then if(DEBUG)print*, $ '** warning ** number of identified '// $ 'couples on plane ',nplx, $ 'exceeds vector dimention ' $ ,'( ',ncouplemax,' )' c good2=.false. c goto 880 !fill ntp and go to next event iflag=1 return endif ncp_plane(nplx) = ncp_plane(nplx) + 1 clx(nplx,ncp_plane(nplx))=icx cly(nply,ncp_plane(nplx))=icy cl_single(icx)=0 cl_single(icy)=0 endif * ---------------------------------------------- 20 continue enddo !end loop on clusters(Y) 10 continue enddo !end loop on clusters(X) do icl=1,nclstr1 if(cl_single(icl).eq.1)then ip=npl(VIEW(icl)) ncls(ip)=ncls(ip)+1 cls(ip,ncls(ip))=icl endif enddo if(DEBUG)then print*,'clusters ',nclstr1 print*,'good ',(cl_good(i),i=1,nclstr1) print*,'singles ',(cl_single(i),i=1,nclstr1) print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes) endif do ip=1,6 ncp_tot=ncp_tot+ncp_plane(ip) enddo c if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!) if(ncp_tot.gt.ncp_max)then if(DEBUG)print*, $ '** warning ** number of identified '// $ 'couples exceeds upper limit for Hough tr. ' $ ,'( ',ncp_max,' )' c good2=.false. c goto 880 !fill ntp and go to next event iflag=1 return endif return end *************************************************** * * * * * * * * * * * * ************************************************** subroutine cl_to_couples_nocharge(iflag) include '../common/commontracker.f' include '../common/common_momanhough.f' include '../common/momanhough_init.f' include '../common/calib.f' include '../common/level1.f' logical DEBUG common/dbg/DEBUG * output flag * -------------- * 0 = good event * 1 = bad event * -------------- integer iflag integer badseed,badcl * init variables ncp_tot=0 do ip=1,nplanes do ico=1,ncouplemax clx(ip,ico)=0 cly(ip,ico)=0 enddo ncp_plane(ip)=0 do icl=1,nclstrmax_level2 cls(ip,icl)=1 enddo ncls(ip)=0 enddo do icl=1,nclstrmax_level2 cl_single(icl)=1 cl_good(icl)=0 enddo * start association ncouples=0 do icx=1,nclstr1 !loop on cluster (X) if(mod(VIEW(icx),2).eq.1)goto 10 * ---------------------------------------------------- * cut on charge (X VIEW) if(dedx(icx).lt.dedx_x_min)then cl_single(icx)=0 goto 10 endif * cut BAD (X VIEW) badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx))) ifirst=INDSTART(icx) if(icx.ne.nclstr1) then ilast=INDSTART(icx+1)-1 else ilast=TOTCLLENGTH endif badcl=badseed do igood=-ngoodstr,ngoodstr ibad=1 if((INDMAX(icx)+igood).gt.ifirst.and. $ (INDMAX(icx)+igood).lt.ilast.and. $ .true.)then ibad=BAD(VIEW(icx), $ nvk(MAXS(icx)+igood), $ nst(MAXS(icx)+igood)) endif badcl=badcl*ibad enddo if(badcl.eq.0)then !<<<<<<<<<<<<<< BAD cut cl_single(icx)=0 !<<<<<<<<<<<<<< BAD cut goto 10 !<<<<<<<<<<<<<< BAD cut endif !<<<<<<<<<<<<<< BAD cut * ---------------------------------------------------- cl_good(icx)=1 nplx=npl(VIEW(icx)) nldx=nld(MAXS(icx),VIEW(icx)) do icy=1,nclstr1 !loop on cluster (Y) if(mod(VIEW(icy),2).eq.0)goto 20 * ---------------------------------------------------- * cut on charge (Y VIEW) if(dedx(icy).lt.dedx_y_min)then cl_single(icy)=0 goto 20 endif * cut BAD (Y VIEW) badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy))) ifirst=INDSTART(icy) if(icy.ne.nclstr1) then ilast=INDSTART(icy+1)-1 else ilast=TOTCLLENGTH endif badcl=badseed do igood=-ngoodstr,ngoodstr ibad=1 if((INDMAX(icy)+igood).gt.ifirst.and. $ (INDMAX(icy)+igood).lt.ilast.and. $ .true.) $ ibad=BAD(VIEW(icy), $ nvk(MAXS(icy)+igood), $ nst(MAXS(icy)+igood)) badcl=badcl*ibad enddo if(badcl.eq.0)then !<<<<<<<<<<<<<< BAD cut cl_single(icy)=0 !<<<<<<<<<<<<<< BAD cut goto 20 !<<<<<<<<<<<<<< BAD cut endif !<<<<<<<<<<<<<< BAD cut * ---------------------------------------------------- cl_good(icy)=1 nply=npl(VIEW(icy)) nldy=nld(MAXS(icy),VIEW(icy)) * ---------------------------------------------- * CONDITION TO FORM A COUPLE * ---------------------------------------------- * geometrical consistency (same plane and ladder) if(nply.eq.nplx.and.nldy.eq.nldx)then * charge correlation * =========================================================== * this version of the subroutine is used for the calibration * thus charge-correlation selection is obviously removed * =========================================================== c$$$ ddd=(dedx(icy) c$$$ $ -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx)) c$$$ ddd=ddd/sqrt(kch(nplx,nldx)**2+1) c$$$ cut=chcut*sch(nplx,nldx) c$$$ if(abs(ddd).gt.cut)goto 20 !charge not consistent * =========================================================== * ------------------> COUPLE <------------------ * check to do not overflow vector dimentions if(ncp_plane(nplx).gt.ncouplemax)then if(DEBUG)print*, $ ' ** warning ** number of identified'// $ ' couples on plane ',nplx, $ ' exceeds vector dimention'// $ ' ( ',ncouplemax,' )' c good2=.false. c goto 880 !fill ntp and go to next event iflag=1 return endif if(ncp_plane(nplx).eq.ncouplemax)then if(DEBUG)print*, $ '** warning ** number of identified '// $ 'couples on plane ',nplx, $ 'exceeds vector dimention ' $ ,'( ',ncouplemax,' )' c good2=.false. c goto 880 !fill ntp and go to next event iflag=1 return endif ncp_plane(nplx) = ncp_plane(nplx) + 1 clx(nplx,ncp_plane(nplx))=icx cly(nply,ncp_plane(nplx))=icy cl_single(icx)=0 cl_single(icy)=0 endif * ---------------------------------------------- 20 continue enddo !end loop on clusters(Y) 10 continue enddo !end loop on clusters(X) do icl=1,nclstr1 if(cl_single(icl).eq.1)then ip=npl(VIEW(icl)) ncls(ip)=ncls(ip)+1 cls(ip,ncls(ip))=icl endif enddo if(DEBUG)then print*,'clusters ',nclstr1 print*,'good ',(cl_good(i),i=1,nclstr1) print*,'singles ',(cl_single(i),i=1,nclstr1) print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes) endif do ip=1,6 ncp_tot=ncp_tot+ncp_plane(ip) enddo c if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!) if(ncp_tot.gt.ncp_max)then if(DEBUG)print*, $ '** warning ** number of identified '// $ 'couples exceeds upper limit for Hough tr. ' $ ,'( ',ncp_max,' )' c good2=.false. c goto 880 !fill ntp and go to next event iflag=1 return endif return end c$$$ subroutine cl_to_couples_2(iflag) c$$$ c$$$ include '../common/commontracker.f' c$$$ include '../common/common_momanhough.f' c$$$ include '../common/momanhough_init.f' c$$$ include '../common/calib.f' c$$$ include '../common/level1.f' c$$$ c$$$ logical DEBUG c$$$ common/dbg/DEBUG c$$$ c$$$* output flag c$$$* -------------- c$$$* 0 = good event c$$$* 1 = bad event c$$$* -------------- c$$$ integer iflag c$$$ c$$$ integer badseed,badcl c$$$ c$$$* init variables c$$$ ncp_tot=0 c$$$ do ip=1,nplanes c$$$ do ico=1,ncouplemax c$$$ clx(ip,ico)=0 c$$$ cly(ip,ico)=0 c$$$ enddo c$$$ ncp_plane(ip)=0 c$$$ do icl=1,nclstrmax_level2 c$$$ cls(ip,icl)=1 c$$$ enddo c$$$ ncls(ip)=0 c$$$ enddo c$$$ do icl=1,nclstrmax_level2 c$$$ cl_single(icl)=1 c$$$ cl_good(icl)=0 c$$$ enddo c$$$ c$$$* start association c$$$ ncouples=0 c$$$ do icx=1,nclstr1 !loop on cluster (X) c$$$ if(mod(VIEW(icx),2).eq.1)goto 10 c$$$ c$$$* ---------------------------------------------------- c$$$* cut on charge (X VIEW) c$$$ if(dedx(icx).lt.dedx_x_min)then c$$$ cl_single(icx)=0 c$$$ goto 10 c$$$ endif c$$$* cut BAD (X VIEW) c$$$ badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx))) c$$$ ifirst=INDSTART(icx) c$$$ if(icx.ne.nclstr1) then c$$$ ilast=INDSTART(icx+1)-1 c$$$ else c$$$ ilast=TOTCLLENGTH c$$$ endif c$$$ badcl=badseed c$$$ do igood=-ngoodstr,ngoodstr c$$$ ibad=1 c$$$ if((INDMAX(icx)+igood).gt.ifirst.and. c$$$ $ (INDMAX(icx)+igood).lt.ilast.and. c$$$ $ .true.)then c$$$ ibad=BAD(VIEW(icx), c$$$ $ nvk(MAXS(icx)+igood), c$$$ $ nst(MAXS(icx)+igood)) c$$$ endif c$$$ badcl=badcl*ibad c$$$ enddo c$$$* print*,'icx ',icx,badcl c$$$ if(badcl.eq.0)then c$$$ cl_single(icx)=0 c$$$ goto 10 c$$$ endif c$$$* ---------------------------------------------------- c$$$ c$$$ cl_good(icx)=1 c$$$ nplx=npl(VIEW(icx)) c$$$ nldx=nld(MAXS(icx),VIEW(icx)) c$$$ c$$$ do icy=1,nclstr1 !loop on cluster (Y) c$$$ if(mod(VIEW(icy),2).eq.0)goto 20 c$$$ c$$$* ---------------------------------------------------- c$$$* cut on charge (Y VIEW) c$$$ if(dedx(icy).lt.dedx_y_min)then c$$$ cl_single(icy)=0 c$$$ goto 20 c$$$ endif c$$$* cut BAD (Y VIEW) c$$$ badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy))) c$$$ ifirst=INDSTART(icy) c$$$ if(icy.ne.nclstr1) then c$$$ ilast=INDSTART(icy+1)-1 c$$$ else c$$$ ilast=TOTCLLENGTH c$$$ endif c$$$ badcl=badseed c$$$ do igood=-ngoodstr,ngoodstr c$$$ ibad=1 c$$$ if((INDMAX(icy)+igood).gt.ifirst.and. c$$$ $ (INDMAX(icy)+igood).lt.ilast.and. c$$$ $ .true.) c$$$ $ ibad=BAD(VIEW(icy), c$$$ $ nvk(MAXS(icy)+igood), c$$$ $ nst(MAXS(icy)+igood)) c$$$ badcl=badcl*ibad c$$$ enddo c$$$* print*,'icy ',icy,badcl c$$$ if(badcl.eq.0)then c$$$ cl_single(icy)=0 c$$$ goto 20 c$$$ endif c$$$* ---------------------------------------------------- c$$$ c$$$ c$$$ cl_good(icy)=1 c$$$ nply=npl(VIEW(icy)) c$$$ nldy=nld(MAXS(icy),VIEW(icy)) c$$$ c$$$* ---------------------------------------------- c$$$* CONDITION TO FORM A COUPLE c$$$* ---------------------------------------------- c$$$* geometrical consistency (same plane and ladder) c$$$ if(nply.eq.nplx.and.nldy.eq.nldx)then c$$$ c$$$c$$$* charge correlation c$$$c$$$ ddd=(dedx(icy) c$$$c$$$ $ -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx)) c$$$c$$$ ddd=ddd/sqrt(kch(nplx,nldx)**2+1) c$$$c$$$ cut=chcut*sch(nplx,nldx) c$$$c$$$ if(abs(ddd).gt.cut)goto 20 !charge not consistent c$$$ c$$$* ------------------> COUPLE <------------------ c$$$* check to do not overflow vector dimentions c$$$ if(ncp_plane(nplx).gt.ncouplemax)then c$$$ if(DEBUG)print*, c$$$ $ ' ** warning ** number of identified'// c$$$ $ ' couples on plane ',nplx, c$$$ $ ' exceeds vector dimention'// c$$$ $ ' ( ',ncouplemax,' )' c$$$c good2=.false. c$$$c goto 880 !fill ntp and go to next event c$$$ iflag=1 c$$$ return c$$$ endif c$$$ c$$$ if(ncp_plane(nplx).eq.ncouplemax)then c$$$ if(DEBUG)print*, c$$$ $ '** warning ** number of identified '// c$$$ $ 'couples on plane ',nplx, c$$$ $ 'exceeds vector dimention ' c$$$ $ ,'( ',ncouplemax,' )' c$$$c good2=.false. c$$$c goto 880 !fill ntp and go to next event c$$$ iflag=1 c$$$ return c$$$ endif c$$$ c$$$ ncp_plane(nplx) = ncp_plane(nplx) + 1 c$$$ clx(nplx,ncp_plane(nplx))=icx c$$$ cly(nply,ncp_plane(nplx))=icy c$$$ cl_single(icx)=0 c$$$ cl_single(icy)=0 c$$$c print*,'couple ',nplx,ncp_plane(nplx),' --- ',icx,icy c$$$ endif c$$$* ---------------------------------------------- c$$$ c$$$ 20 continue c$$$ enddo !end loop on clusters(Y) c$$$ c$$$ 10 continue c$$$ enddo !end loop on clusters(X) c$$$ c$$$ c$$$ do icl=1,nclstr1 c$$$ if(cl_single(icl).eq.1)then c$$$ ip=npl(VIEW(icl)) c$$$ ncls(ip)=ncls(ip)+1 c$$$ cls(ip,ncls(ip))=icl c$$$ endif c$$$ enddo c$$$ c$$$ c$$$ if(DEBUG)then c$$$ print*,'clusters ',nclstr1 c$$$ print*,'good ',(cl_good(i),i=1,nclstr1) c$$$ print*,'singles ',(cl_single(i),i=1,nclstr1) c$$$ print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes) c$$$ endif c$$$ c$$$ do ip=1,6 c$$$ ncp_tot=ncp_tot+ncp_plane(ip) c$$$ enddo c$$$c if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!) c$$$ c$$$ if(ncp_tot.gt.ncp_max)then c$$$ if(DEBUG)print*, c$$$ $ '** warning ** number of identified '// c$$$ $ 'couples exceeds upper limit for Hough tr. ' c$$$ $ ,'( ',ncp_max,' )' c$$$c good2=.false. c$$$c goto 880 !fill ntp and go to next event c$$$ iflag=1 c$$$ return c$$$ endif c$$$ c$$$ return c$$$ end *************************************************** * * * * * * * * * * * * ************************************************** subroutine cp_to_doubtrip(iflag) c***************************************************** c 02/02/2006 modified by Elena Vannuccini --> (1) c***************************************************** include '../common/commontracker.f' include '../common/common_momanhough.f' include '../common/momanhough_init.f' include '../common/common_xyzPAM.f' include '../common/common_mini_2.f' include '../common/calib.f' include '../common/level1.f' logical DEBUG common/dbg/DEBUG * output flag * -------------- * 0 = good event * 1 = bad event * -------------- integer iflag * ----------------------------- * DOUBLETS/TRIPLETS coordinates c double precision xm1,ym1,zm1 c double precision xm2,ym2,zm2 c double precision xm3,ym3,zm3 real xm1,ym1,zm1 real xm2,ym2,zm2 real xm3,ym3,zm3 * ----------------------------- * variable needed for tricircle: real xp(3),zp(3)!TRIPLETS coordinates, to find a circle EQUIVALENCE (xm1,xp(1)) EQUIVALENCE (xm2,xp(2)) EQUIVALENCE (xm3,xp(3)) EQUIVALENCE (zm1,zp(1)) EQUIVALENCE (zm2,zp(2)) EQUIVALENCE (zm3,zp(3)) real angp(3),resp(3),chi real xc,zc,radius * ----------------------------- ndblt=0 !number of doublets ntrpt=0 !number of triplets do ip1=1,(nplanes-1) !loop on planes - COPPIA 1 do is1=1,2 !loop on sensors - COPPIA 1 do icp1=1,ncp_plane(ip1) !loop on COPPIA 1 icx1=clx(ip1,icp1) icy1=cly(ip1,icp1) c call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1) call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1) xm1=xPAM ym1=yPAM zm1=zPAM c print*,'***',is1,xm1,ym1,zm1 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2 do is2=1,2 !loop on sensors -ndblt COPPIA 2 do icp2=1,ncp_plane(ip2) !loop on COPPIA 2 icx2=clx(ip2,icp2) icy2=cly(ip2,icp2) c call xyz_PAM c $ (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1) call xyz_PAM $ (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1) xm2=xPAM ym2=yPAM zm2=zPAM * - - - - - - - - - - - - - - - - - - - - - - - - - - - - * track parameters on Y VIEW * (2 couples needed) * - - - - - - - - - - - - - - - - - - - - - - - - - - - - if(ndblt.eq.ndblt_max)then if(DEBUG)print*, $ '** warning ** number of identified '// $ 'doublets exceeds vector dimention ' $ ,'( ',ndblt_max,' )' c good2=.false. c goto 880 !fill ntp and go to next event iflag=1 return endif ndblt = ndblt + 1 * store doublet info cpyz1(ndblt)=id_cp(ip1,icp1,is1) cpyz2(ndblt)=id_cp(ip2,icp2,is2) * tg(th_yz) alfayz2(ndblt)=(ym1-ym2)/(zm1-zm2) * y0 (cm) alfayz1(ndblt)=alfayz2(ndblt)*(zini-zm1)+ym1 **** -----------------------------------------------**** **** reject non phisical couples **** **** -----------------------------------------------**** if( $ abs(alfayz2(ndblt)).gt.alfyz2_max $ .or. $ abs(alfayz1(ndblt)).gt.alfyz1_max $ )ndblt = ndblt-1 c$$$ if(iev.eq.33)then c$$$ print*,'********* ',ndblt,' -- ',icp1,icp2,is1,is2 c$$$ $ ,' || ',icx1,icy1,icx2,icy2 c$$$ $ ,' || ',xm1,ym1,xm2,ym2 c$$$ $ ,' || ',alfayz2(ndblt),alfayz1(ndblt) c$$$ endif c$$$ * - - - - - - - - - - - - - - - - - - - - - - - - - - - - * track parameters on Y VIEW - end * - - - - - - - - - - - - - - - - - - - - - - - - - - - - if(ip2.eq.nplanes)goto 30 !no possible combination with 3 couples do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3 do is3=1,2 !loop on sensors - COPPIA 3 do icp3=1,ncp_plane(ip3) !loop on COPPIA 3 icx3=clx(ip3,icp3) icy3=cly(ip3,icp3) c call xyz_PAM c $ (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1) call xyz_PAM $ (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1) xm3=xPAM ym3=yPAM zm3=zPAM * find the circle passing through the three points call tricircle(3,xp,zp,angp,resp,chi $ ,xc,zc,radius,iflag) c print*,xc,zc,radius * the circle must intersect the reference plane if( c $ (xc.le.-1.*xclimit.or. c $ xc.ge.xclimit).and. $ radius**2.ge.(ZINI-zc)**2.and. $ iflag.eq.0.and. $ .true.)then * - - - - - - - - - - - - - - - - - - - - - - - - - - - - * track parameters on X VIEW * (3 couples needed) * - - - - - - - - - - - - - - - - - - - - - - - - - - - - if(ntrpt.eq.ntrpt_max)then if(DEBUG)print*, $ '** warning ** number of identified '// $ 'triplets exceeds vector dimention ' $ ,'( ',ntrpt_max,' )' c good2=.false. c goto 880 !fill ntp and go to next event iflag=1 return endif ntrpt = ntrpt +1 * store triplet info cpxz1(ntrpt)=id_cp(ip1,icp1,is1) cpxz2(ntrpt)=id_cp(ip2,icp2,is2) cpxz3(ntrpt)=id_cp(ip3,icp3,is3) if(xc.lt.0)then *************POSITIVE DEFLECTION alfaxz1(ntrpt) = xc+sqrt(radius**2-(ZINI-zc)**2) alfaxz2(ntrpt) = (ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2) alfaxz3(ntrpt) = 1/radius else *************NEGATIVE DEFLECTION alfaxz1(ntrpt) = xc-sqrt(radius**2-(ZINI-zc)**2) alfaxz2(ntrpt) = -(ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2) alfaxz3(ntrpt) = -1/radius endif **** -----------------------------------------------**** **** reject non phisical triplets **** **** -----------------------------------------------**** if( $ abs(alfaxz2(ntrpt)).gt.alfxz2_max $ .or. $ abs(alfaxz1(ntrpt)).gt.alfxz1_max $ )ntrpt = ntrpt-1 c print*,alfaxz1(ntrpt),alfaxz2(ntrpt),alfaxz3(ntrpt) * - - - - - - - - - - - - - - - - - - - - - - - - - - - - * track parameters on X VIEW - end * - - - - - - - - - - - - - - - - - - - - - - - - - - - - endif enddo !end loop on COPPIA 3 enddo !end loop on sensors - COPPIA 3 enddo !end loop on planes - COPPIA 3 30 continue 1 enddo !end loop on COPPIA 2 enddo !end loop on sensors - COPPIA 2 enddo !end loop on planes - COPPIA 2 enddo !end loop on COPPIA1 enddo !end loop on sensors - COPPIA 1 enddo !end loop on planes - COPPIA 1 if(DEBUG)then print*,'--- doublets ',ndblt print*,'--- triplets ',ntrpt endif c goto 880 !ntp fill return end *************************************************** * * * * * * * * * * * * ************************************************** subroutine doub_to_YZcloud(iflag) include '../common/commontracker.f' include '../common/common_momanhough.f' include '../common/momanhough_init.f' logical DEBUG common/dbg/DEBUG * output flag * -------------- * 0 = good event * 1 = bad event * -------------- integer iflag integer db_used(ndblt_max) integer db_temp(ndblt_max) integer db_all(ndblt_max) !stores db ID in each cloud integer hit_plane(nplanes) * mask for used couples integer cp_useds1(ncouplemaxtot) ! sensor 1 integer cp_useds2(ncouplemaxtot) ! sensor 2 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ * classification of DOUBLETS * according to distance in parameter space * (cloud = group of points (doublets) in parameter space) *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ do idb=1,ndblt db_used(idb)=0 enddo distance=0 nclouds_yz=0 !number of clouds npt_tot=0 do idb1=1,ndblt !loop (1) on DOUBLETS if(db_used(idb1).eq.1)goto 2228 !db already included in a cloud c print*,'--------------' c print*,'** ',idb1,' **' do icp=1,ncp_tot cp_useds1(icp)=0 !init cp_useds2(icp)=0 !init enddo do idb=1,ndblt db_all(idb)=0 enddo if(cpyz1(idb1).gt.0)cp_useds2(cpyz1(idb1))=1 if(cpyz1(idb1).lt.0)cp_useds1(-cpyz1(idb1))=1 if(cpyz2(idb1).gt.0)cp_useds2(cpyz2(idb1))=1 if(cpyz2(idb1).lt.0)cp_useds1(-cpyz2(idb1))=1 temp1 = alfayz1(idb1) temp2 = alfayz2(idb1) npt=1 !counter of points in the cloud db_all(npt) = idb1 nptloop=1 db_temp(1)=idb1 88 continue npv=0 !# new points inlcuded do iloop=1,nptloop idbref=db_temp(iloop) !local point of reference ccccc if(db_used(idbref).eq.1)goto 1188 !next do idb2=1,ndblt !loop (2) on DOUBLETS if(idb2.eq.idbref)goto 1118 !next doublet if(db_used(idb2).eq.1)goto 1118 * doublet distance in parameter space distance= $ ((alfayz1(idbref)-alfayz1(idb2))/Dalfayz1)**2 $ +((alfayz2(idbref)-alfayz2(idb2))/Dalfayz2)**2 distance = sqrt(distance) c$$$ if(iev.eq.33)then c$$$ if(distance.lt.100) c$$$ $ print*,'********* ',idb1,idbref,idb2,distance c$$$ if(distance.lt.100) c$$$ $ print*,'********* ',alfayz1(idbref),alfayz1(idb2) c$$$ $ ,alfayz2(idbref),alfayz2(idb2) c$$$ endif if(distance.lt.cutdistyz)then c print*,idb1,idb2,distance,' cloud ',nclouds_yz if(cpyz1(idb2).gt.0)cp_useds2(cpyz1(idb2))=1 if(cpyz1(idb2).lt.0)cp_useds1(-cpyz1(idb2))=1 if(cpyz2(idb2).gt.0)cp_useds2(cpyz2(idb2))=1 if(cpyz2(idb2).lt.0)cp_useds1(-cpyz2(idb2))=1 npt = npt + 1 !counter of points in the cloud npv = npv +1 db_temp(npv) = idb2 db_used(idbref) = 1 db_used(idb2) = 1 db_all(npt) = idb2 temp1 = temp1 + alfayz1(idb2) temp2 = temp2 + alfayz2(idb2) c print*,'* idbref,idb2 ',idbref,idb2 endif 1118 continue enddo !end loop (2) on DOUBLETS 1188 continue enddo !end loop on... bo? nptloop=npv if(nptloop.ne.0)goto 88 * ------------------------------------------ * stores the cloud only if * 1) it includes a minimum number of REAL couples * 1bis) it inlcudes a minimum number of doublets * 2) it is not already stored * ------------------------------------------ do ip=1,nplanes hit_plane(ip)=0 enddo ncpused=0 do icp=1,ncp_tot if(cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0)then ncpused=ncpused+1 ip=ip_cp(icp) hit_plane(ip)=1 endif enddo nplused=0 do ip=1,nplanes nplused=nplused+ hit_plane(ip) enddo c print*,'>>>> ',ncpused,npt,nplused if(ncpused.lt.ncpyz_min)goto 2228 !next doublet if(npt.lt.nptyz_min)goto 2228 !next doublet if(nplused.lt.nplyz_min)goto 2228 !next doublet * ~~~~~~~~~~~~~~~~~ * >>> NEW CLOUD <<< if(nclouds_yz.ge.ncloyz_max)then if(DEBUG)print*, $ '** warning ** number of identified '// $ 'YZ clouds exceeds vector dimention ' $ ,'( ',ncloyz_max,' )' c good2=.false. c goto 880 !fill ntp and go to next event iflag=1 return endif nclouds_yz = nclouds_yz + 1 !increase counter alfayz1_av(nclouds_yz) = temp1/npt !store average parameter alfayz2_av(nclouds_yz) = temp2/npt ! " do icp=1,ncp_tot cpcloud_yz(nclouds_yz,icp)= $ cp_useds1(icp)+2*cp_useds2(icp) !store cp info enddo ptcloud_yz(nclouds_yz)=npt c ptcloud_yz_nt(nclouds_yz)=npt do ipt=1,npt db_cloud(npt_tot+ipt) = db_all(ipt) c print*,'>> ',ipt,db_all(ipt) enddo npt_tot=npt_tot+npt if(DEBUG)then print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~' print*,'>>>> cloud ',nclouds_yz,' --- ',npt,' points' print*,'- alfayz1 ',alfayz1_av(nclouds_yz) print*,'- alfayz2 ',alfayz2_av(nclouds_yz) print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot) print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot) print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes) c$$$ print*,'nt-uple: ptcloud_yz(',nclouds_yz,') = ' c$$$ $ ,ptcloud_yz(nclouds_yz) c$$$ print*,'nt-uple: db_cloud(...) = ' c$$$ $ ,(db_cloud(iii),iii=npt_tot-npt+1,npt_tot) endif * >>> NEW CLOUD <<< * ~~~~~~~~~~~~~~~~~ 2228 continue enddo !end loop (1) on DOUBLETS if(DEBUG)then print*,'---------------------- ' print*,'Y-Z total clouds ',nclouds_yz print*,' ' endif return end *************************************************** * * * * * * * * * * * * ************************************************** subroutine trip_to_XZcloud(iflag) include '../common/commontracker.f' include '../common/common_momanhough.f' include '../common/momanhough_init.f' logical DEBUG common/dbg/DEBUG * output flag * -------------- * 0 = good event * 1 = bad event * -------------- integer iflag integer tr_used(ntrpt_max) integer tr_temp(ntrpt_max) integer tr_incl(ntrpt_max) integer tr_all(ntrpt_max) !stores tr ID in each cloud integer hit_plane(nplanes) * mask for used couples integer cp_useds1(ncouplemaxtot) ! sensor 1 integer cp_useds2(ncouplemaxtot) ! sensor 2 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ * classification of TRIPLETS * according to distance in parameter space *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ do itr=1,ntrpt tr_used(itr)=0 enddo distance=0 nclouds_xz=0 !number of clouds npt_tot=0 !total number of selected triplets do itr1=1,ntrpt !loop (1) on TRIPLETS if(tr_used(itr1).eq.1)goto 22288 !already included in a cloud c print*,'--------------' c print*,'** ',itr1,' **' do icp=1,ncp_tot cp_useds1(icp)=0 cp_useds2(icp)=0 enddo do itr=1,ntrpt tr_all(itr)=0 !list of included triplets enddo if(cpxz1(itr1).gt.0)cp_useds2(cpxz1(itr1))=1 if(cpxz1(itr1).lt.0)cp_useds1(-cpxz1(itr1))=1 if(cpxz2(itr1).gt.0)cp_useds2(cpxz2(itr1))=1 if(cpxz2(itr1).lt.0)cp_useds1(-cpxz2(itr1))=1 if(cpxz3(itr1).gt.0)cp_useds2(cpxz3(itr1))=1 if(cpxz3(itr1).lt.0)cp_useds1(-cpxz3(itr1))=1 temp1 = alfaxz1(itr1) temp2 = alfaxz2(itr1) temp3 = alfaxz3(itr1) npt=1 !counter of points in the cloud tr_all(npt) = itr1 nptloop=1 c tr_temp(1)=itr1 tr_incl(1)=itr1 8881 continue npv=0 !# new points inlcuded do iloop=1,nptloop itrref=tr_incl(iloop) !local point of reference do itr2=1,ntrpt !loop (2) on TRIPLETS if(itr2.eq.itr1)goto 11188 !next triplet if(tr_used(itr2).eq.1)goto 11188 !next triplet * triplet distance in parameter space * solo i due parametri spaziali per il momemnto distance= $ ((alfaxz1(itrref)-alfaxz1(itr2))/Dalfaxz1)**2 $ +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2 distance = sqrt(distance) if(distance.lt.cutdistxz)then c print*,idb1,idb2,distance,' cloud ',nclouds_yz if(cpxz1(itr2).gt.0)cp_useds2(cpxz1(itr2))=1 if(cpxz1(itr2).lt.0)cp_useds1(-cpxz1(itr2))=1 if(cpxz2(itr2).gt.0)cp_useds2(cpxz2(itr2))=1 if(cpxz2(itr2).lt.0)cp_useds1(-cpxz2(itr2))=1 if(cpxz3(itr2).gt.0)cp_useds2(cpxz3(itr2))=1 if(cpxz3(itr2).lt.0)cp_useds1(-cpxz3(itr2))=1 npt = npt + 1 !counter of points in the cloud npv = npv +1 tr_temp(npv) = itr2 tr_used(itrref) = 1 tr_used(itr2) = 1 tr_all(npt) = itr2 temp1 = temp1 + alfaxz1(itr2) temp2 = temp2 + alfaxz2(itr2) temp3 = temp3 + alfaxz3(itr2) c print*,'* itrref,itr2 ',itrref,itr2,distance endif 11188 continue enddo !end loop (2) on TRIPLETS 11888 continue enddo !end loop on... bo? nptloop=npv do i=1,npv tr_incl(i)=tr_temp(i) enddo if(nptloop.ne.0)goto 8881 * ------------------------------------------ * stores the cloud only if * 1) it includes a minimum number of REAL couples * 1bis) * 2) it is not already stored * ------------------------------------------ c print*,'check cp_used' do ip=1,nplanes hit_plane(ip)=0 enddo ncpused=0 do icp=1,ncp_tot if(cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0)then ncpused=ncpused+1 ip=ip_cp(icp) hit_plane(ip)=1 endif enddo nplused=0 do ip=1,nplanes nplused=nplused+ hit_plane(ip) enddo if(ncpused.lt.ncpxz_min)goto 22288 !next triplet if(npt.lt.nptxz_min)goto 22288 !next triplet if(nplused.lt.nplxz_min)goto 22288 !next doublet * ~~~~~~~~~~~~~~~~~ * >>> NEW CLOUD <<< if(nclouds_xz.ge.ncloxz_max)then if(DEBUG)print*, $ '** warning ** number of identified '// $ 'XZ clouds exceeds vector dimention ' $ ,'( ',ncloxz_max,' )' c good2=.false. c goto 880 !fill ntp and go to next event iflag=1 return endif nclouds_xz = nclouds_xz + 1 !increase counter alfaxz1_av(nclouds_xz) = temp1/npt !store average parameter alfaxz2_av(nclouds_xz) = temp2/npt ! " alfaxz3_av(nclouds_xz) = temp3/npt ! " do icp=1,ncp_tot cpcloud_xz(nclouds_xz,icp)= $ cp_useds1(icp)+2*cp_useds2(icp) !store cp info enddo ptcloud_xz(nclouds_xz)=npt do ipt=1,npt tr_cloud(npt_tot+ipt) = tr_all(ipt) enddo npt_tot=npt_tot+npt if(DEBUG)then print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~' print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points' print*,'- alfaxz1 ',alfaxz1_av(nclouds_xz) print*,'- alfaxz2 ',alfaxz2_av(nclouds_xz) print*,'- alfaxz3 ',alfaxz3_av(nclouds_xz) print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot) print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot) print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes) c$$$ print*,'nt-uple: ptcloud_xz(',nclouds_xz,') = ' c$$$ $ ,ptcloud_xz(nclouds_xz) c$$$ print*,'nt-uple: tr_cloud(...) = ' c$$$ $ ,(tr_cloud(iii),iii=npt_tot-npt+1,npt_tot) endif * >>> NEW CLOUD <<< * ~~~~~~~~~~~~~~~~~ 22288 continue enddo !end loop (1) on DOUBLETS if(DEBUG)then print*,'---------------------- ' print*,'X-Z total clouds ',nclouds_xz print*,' ' endif return end *************************************************** * * * * * * * * * * * * ************************************************** subroutine clouds_to_ctrack(iflag) c***************************************************** c 02/02/2006 modified by Elena Vannuccini --> (1) c***************************************************** include '../common/commontracker.f' include '../common/common_momanhough.f' include '../common/common_xyzPAM.f' include '../common/common_mini_2.f' include '../common/common_mech.f' include '../common/momanhough_init.f' logical DEBUG common/dbg/DEBUG * output flag * -------------- * 0 = good event * 1 = bad event * -------------- integer iflag * ----------------------------------------------------------- * mask to store (locally) the couples included * in the intersection bewteen a XZ and YZ cloud integer cpintersec(ncouplemaxtot) * ----------------------------------------------------------- * list of matching couples in the combination * between a XZ and YZ cloud integer cp_match(nplanes,ncouplemax) integer ncp_match(nplanes) * ----------------------------------------------------------- integer hit_plane(nplanes) * ----------------------------------------------------------- * variables for track fitting double precision AL_INI(5) double precision tath * ----------------------------------------------------------- c real fitz(nplanes) !z coordinates of the planes in cm ntracks=0 !counter of track candidates do iyz=1,nclouds_yz !loop on YZ couds do ixz=1,nclouds_xz !loop on XZ couds * -------------------------------------------------- * check of consistency of the clouds * ---> required a minimum number of matching couples * the track fit will be performed on the INTERSECTION * of the two clouds * -------------------------------------------------- do ip=1,nplanes hit_plane(ip)=0 ncp_match(ip)=0 do icpp=1,ncouplemax cp_match(ip,icpp)=0 !init couple list enddo enddo ncp_ok=0 do icp=1,ncp_tot !loop on couples * get info on cpintersec(icp)=min( $ cpcloud_yz(iyz,icp), $ cpcloud_xz(ixz,icp)) if( $ (cpcloud_yz(iyz,icp).eq.1.and.cpcloud_xz(ixz,icp).eq.2).or. $ (cpcloud_yz(iyz,icp).eq.2.and.cpcloud_xz(ixz,icp).eq.1).or. $ .false.)cpintersec(icp)=0 if(cpintersec(icp).ne.0)then ncp_ok=ncp_ok+1 ip=ip_cp(icp) hit_plane(ip)=1 if(cpintersec(icp).eq.1)then * 1) only the couple image in sensor 1 matches id=-icp ncp_match(ip)=ncp_match(ip)+1 cp_match(ip,ncp_match(ip))=id elseif(cpintersec(icp).eq.2)then * 2) only the couple image in sensor 2 matches id=icp ncp_match(ip)=ncp_match(ip)+1 cp_match(ip,ncp_match(ip))=id else * 3) both couple images match id=icp do is=1,2 id=-id ncp_match(ip)=ncp_match(ip)+1 cp_match(ip,ncp_match(ip))=id enddo endif endif !end matching condition enddo !end loop on couples nplused=0 do ip=1,nplanes nplused=nplused+ hit_plane(ip) enddo if(nplused.lt.nplxz_min)goto 888 !next doublet if(ncp_ok.lt.ncpok)goto 888 !next cloud if(DEBUG)then print*,'Combination ',iyz,ixz $ ,' db ',ptcloud_yz(iyz) $ ,' tr ',ptcloud_xz(ixz) $ ,' -----> # matching couples ',ncp_ok endif c$$$ print*,'~~~~~~~~~~~~~~~~~~~~~~~~~' c$$$ print*,'Configurazione cluster XZ' c$$$ print*,'1 -- ',(clx(1,i),i=1,ncp_plane(1)) c$$$ print*,'2 -- ',(clx(2,i),i=1,ncp_plane(1)) c$$$ print*,'3 -- ',(clx(3,i),i=1,ncp_plane(1)) c$$$ print*,'4 -- ',(clx(4,i),i=1,ncp_plane(1)) c$$$ print*,'5 -- ',(clx(5,i),i=1,ncp_plane(1)) c$$$ print*,'6 -- ',(clx(6,i),i=1,ncp_plane(1)) c$$$ print*,'Configurazione cluster YZ' c$$$ print*,'1 -- ',(cly(1,i),i=1,ncp_plane(1)) c$$$ print*,'2 -- ',(cly(2,i),i=1,ncp_plane(1)) c$$$ print*,'3 -- ',(cly(3,i),i=1,ncp_plane(1)) c$$$ print*,'4 -- ',(cly(4,i),i=1,ncp_plane(1)) c$$$ print*,'5 -- ',(cly(5,i),i=1,ncp_plane(1)) c$$$ print*,'6 -- ',(cly(6,i),i=1,ncp_plane(1)) c$$$ print*,'~~~~~~~~~~~~~~~~~~~~~~~~~' * -------> INITIAL GUESS <------- AL_INI(1)=dreal(alfaxz1_av(ixz)) AL_INI(2)=dreal(alfayz1_av(iyz)) AL_INI(4)=datan(dreal(alfayz2_av(iyz)) $ /dreal(alfaxz2_av(ixz))) tath=-dreal(alfaxz2_av(ixz))/dcos(AL_INI(4)) AL_INI(3)=tath/sqrt(1+tath**2) AL_INI(5)=(1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0. c print*,'*******',AL_INI(5) if(AL_INI(5).gt.defmax)goto 888 !next cloud c print*,'alfaxz2, alfayz2 ' c $ ,alfaxz2_av(ixz),alfayz2_av(iyz) * -------> INITIAL GUESS <------- c print*,'AL_INI ',(al_ini(i),i=1,5) if(DEBUG)then print*,'1 >>> ',(cp_match(6,i),i=1,ncp_match(6)) print*,'2 >>> ',(cp_match(5,i),i=1,ncp_match(5)) print*,'3 >>> ',(cp_match(4,i),i=1,ncp_match(4)) print*,'4 >>> ',(cp_match(3,i),i=1,ncp_match(3)) print*,'5 >>> ',(cp_match(2,i),i=1,ncp_match(2)) print*,'6 >>> ',(cp_match(1,i),i=1,ncp_match(1)) endif do icp1=1,max(1,ncp_match(1)) hit_plane(1)=icp1 if(ncp_match(1).eq.0)hit_plane(1)=0 !-icp1 do icp2=1,max(1,ncp_match(2)) hit_plane(2)=icp2 if(ncp_match(2).eq.0)hit_plane(2)=0 !-icp2 do icp3=1,max(1,ncp_match(3)) hit_plane(3)=icp3 if(ncp_match(3).eq.0)hit_plane(3)=0 !-icp3 do icp4=1,max(1,ncp_match(4)) hit_plane(4)=icp4 if(ncp_match(4).eq.0)hit_plane(4)=0 !-icp4 do icp5=1,max(1,ncp_match(5)) hit_plane(5)=icp5 if(ncp_match(5).eq.0)hit_plane(5)=0 !-icp5 do icp6=1,max(1,ncp_match(6)) hit_plane(6)=icp6 if(ncp_match(6).eq.0)hit_plane(6)=0 !-icp6 call track_init !init TRACK common do ip=1,nplanes !loop on planes if(hit_plane(ip).ne.0)then id=cp_match(ip,hit_plane(ip)) is=is_cp(id) icp=icp_cp(id) if(ip_cp(id).ne.ip) $ print*,'OKKIO!!' $ ,'id ',id,is,icp $ ,ip_cp(id),ip icx=clx(ip,icp) icy=cly(ip,icp) * ************************* c call xyz_PAM(icx,icy,is, c $ 'COG2','COG2',0.,0.) call xyz_PAM(icx,icy,is, !(1) $ PFAdef,PFAdef,0.,0.) !(1) * ************************* * ----------------------------- xgood(nplanes-ip+1)=1. ygood(nplanes-ip+1)=1. xm(nplanes-ip+1)=xPAM ym(nplanes-ip+1)=yPAM zm(nplanes-ip+1)=zPAM resx(nplanes-ip+1)=resxPAM resy(nplanes-ip+1)=resyPAM * ----------------------------- endif enddo !end loop on planes * ********************************************************** * ************************** FIT *** FIT *** FIT *** FIT *** * ********************************************************** do i=1,5 AL(i)=AL_INI(i) enddo ifail=0 !error flag in chi^2 computation jstep=0 !number of minimization steps call mini_2(jstep,ifail) if(ifail.ne.0) then if(DEBUG)then print *, $ '*** MINIMIZATION FAILURE *** ' $ //'(mini_2 in clouds_to_ctrack)' endif chi2=-chi2 endif * ********************************************************** * ************************** FIT *** FIT *** FIT *** FIT *** * ********************************************************** if(chi2.le.0.)goto 666 * -------------------------- * STORE candidate TRACK INFO * -------------------------- if(ntracks.eq.NTRACKSMAX)then if(DEBUG)print*, $ '** warning ** number of candidate tracks '// $ ' exceeds vector dimension ' $ ,'( ',NTRACKSMAX,' )' c good2=.false. c goto 880 !fill ntp and go to next event iflag=1 return endif ntracks = ntracks + 1 c$$$ ndof=0 do ip=1,nplanes c$$$ ndof=ndof c$$$ $ +int(xgood(ip)) c$$$ $ +int(ygood(ip)) XV_STORE(ip,ntracks)=sngl(xv(ip)) YV_STORE(ip,ntracks)=sngl(yv(ip)) ZV_STORE(ip,ntracks)=sngl(zv(ip)) XM_STORE(ip,ntracks)=sngl(xm(ip)) YM_STORE(ip,ntracks)=sngl(ym(ip)) ZM_STORE(ip,ntracks)=sngl(zm(ip)) RESX_STORE(ip,ntracks)=sngl(resx(ip)) RESY_STORE(ip,ntracks)=sngl(resy(ip)) XV_STORE(ip,ntracks)=sngl(xv(ip)) YV_STORE(ip,ntracks)=sngl(yv(ip)) ZV_STORE(ip,ntracks)=sngl(zv(ip)) AXV_STORE(ip,ntracks)=sngl(axv(ip)) AYV_STORE(ip,ntracks)=sngl(ayv(ip)) XGOOD_STORE(ip,ntracks)=sngl(xgood(ip)) YGOOD_STORE(ip,ntracks)=sngl(ygood(ip)) if(hit_plane(ip).ne.0)then CP_STORE(nplanes-ip+1,ntracks)= $ cp_match(ip,hit_plane(ip)) else CP_STORE(nplanes-ip+1,ntracks)=0 endif CLS_STORE(nplanes-ip+1,ntracks)=0 do i=1,5 AL_STORE(i,ntracks)=sngl(AL(i)) enddo enddo c$$$ * Number of Degree Of Freedom c$$$ ndof=ndof-5 c$$$ * reduced chi^2 c$$$ rchi2=chi2/dble(ndof) RCHI2_STORE(ntracks)=chi2 * -------------------------------- * STORE candidate TRACK INFO - end * -------------------------------- 666 continue enddo !end loop on cp in plane 6 enddo !end loop on cp in plane 5 enddo !end loop on cp in plane 4 enddo !end loop on cp in plane 3 enddo !end loop on cp in plane 2 enddo !end loop on cp in plane 1 888 continue enddo !end loop on XZ couds enddo !end loop on YZ couds if(ntracks.eq.0)then iflag=1 return endif if(DEBUG)then print*,'****** TRACK CANDIDATES ***********' print*,'# R. chi2 RIG' do i=1,ntracks print*,i,' --- ',rchi2_store(i),' --- ' $ ,1./abs(AL_STORE(5,i)) enddo print*,'***********************************' endif return end *************************************************** * * * * * * * * * * * * ************************************************** subroutine refine_track(ibest) c****************************************************** cccccc 06/10/2005 modified by elena vannuccini ---> (1) cccccc 31/01/2006 modified by elena vannuccini ---> (2) c****************************************************** include '../common/commontracker.f' include '../common/common_momanhough.f' include '../common/common_xyzPAM.f' include '../common/common_mini_2.f' include '../common/common_mech.f' include '../common/momanhough_init.f' include '../common/level1.f' include '../common/calib.f' logical DEBUG common/dbg/DEBUG * flag to chose PFA character*10 PFA common/FINALPFA/PFA * ================================================= * new estimate of positions using ETA algorithm * and * search for new couples and single clusters to add * ================================================= call track_init do ip=1,nplanes !loop on planes * ------------------------------------------------- * If the plane has been already included, it just * computes again the coordinates of the x-y couple * using improved PFAs * ------------------------------------------------- if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and. $ YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then id=CP_STORE(nplanes-ip+1,ibest) is=is_cp(id) icp=icp_cp(id) if(ip_cp(id).ne.ip) $ print*,'OKKIO!!' $ ,'id ',id,is,icp $ ,ip_cp(id),ip icx=clx(ip,icp) icy=cly(ip,icp) call xyz_PAM(icx,icy,is, c $ 'ETA2','ETA2', $ PFA,PFA, $ AXV_STORE(nplanes-ip+1,ibest), $ AYV_STORE(nplanes-ip+1,ibest)) c$$$ call xyz_PAM(icx,icy,is, c$$$ $ 'COG2','COG2', c$$$ $ 0., c$$$ $ 0.) xm(nplanes-ip+1) = xPAM ym(nplanes-ip+1) = yPAM zm(nplanes-ip+1) = zPAM xgood(nplanes-ip+1) = 1 ygood(nplanes-ip+1) = 1 resx(nplanes-ip+1) = resxPAM resy(nplanes-ip+1) = resyPAM c dedxtrk(nplanes-ip+1) = (dedx(icx)+dedx(icy))/2. !(1) dedxtrk_x(nplanes-ip+1)=dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2) dedxtrk_y(nplanes-ip+1)=dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2) * ------------------------------------------------- * If the plane has NOT been already included, * it tries to include a COUPLE or a single cluster * ------------------------------------------------- else xgood(nplanes-ip+1)=0 ygood(nplanes-ip+1)=0 * -------------------------------------------------------------- * determine which ladder and sensor are intersected by the track xP=XV_STORE(nplanes-ip+1,ibest) yP=YV_STORE(nplanes-ip+1,ibest) zP=ZV_STORE(nplanes-ip+1,ibest) call whichsensor(ip,xP,yP,nldt,ist) * if the track hit the plane in a dead area, go to the next plane if(nldt.eq.0.or.ist.eq.0)goto 133 * -------------------------------------------------------------- if(DEBUG)then print*, $ '------ Plane ',ip,' intersected on LADDER ',nldt $ ,' SENSOR ',ist print*, $ '------ coord: ',XP,YP endif * =========================================== * STEP 1 >>>>>>> try to include a new couple * =========================================== c if(DEBUG)print*,'>>>> try to include a new couple' distmin=1000000. xmm = 0. ymm = 0. zmm = 0. rxmm = 0. rymm = 0. dedxmmx = 0. !(1) dedxmmy = 0. !(1) idm = 0 !ID of the closer couple distance=0. do icp=1,ncp_plane(ip) !loop on couples on plane icp icx=clx(ip,icp) icy=cly(ip,icp) if(LADDER(icx).ne.nldt.or. !If the ladder number does not match $ cl_used(icx).eq.1.or. !or the X cluster is already used $ cl_used(icy).eq.1.or. !or the Y cluster is already used $ .false.)goto 1188 !then jump to next couple. * call xyz_PAM(icx,icy,ist, $ PFA,PFA, c $ 'ETA2','ETA2', $ AXV_STORE(nplanes-ip+1,ibest), $ AYV_STORE(nplanes-ip+1,ibest)) distance = distance_to(XP,YP) id=id_cp(ip,icp,ist) if(DEBUG)print*,'( couple ',id $ ,' ) normalized distance ',distance if(distance.lt.distmin)then xmm = xPAM ymm = yPAM zmm = zPAM rxmm = resxPAM rymm = resyPAM distmin = distance idm = id c dedxmm = (dedx(icx)+dedx(icy))/2. !(1) dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2) dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2) endif 1188 continue enddo !end loop on couples on plane icp if(distmin.le.clinc)then * ----------------------------------- xm(nplanes-ip+1) = xmm !<<< ym(nplanes-ip+1) = ymm !<<< zm(nplanes-ip+1) = zmm !<<< xgood(nplanes-ip+1) = 1 !<<< ygood(nplanes-ip+1) = 1 !<<< resx(nplanes-ip+1)=rxmm !<<< resy(nplanes-ip+1)=rymm !<<< c dedxtrk(nplanes-ip+1) = dedxmm !<<< !(1) dedxtrk_x(nplanes-ip+1) = dedxmmx !(1) dedxtrk_y(nplanes-ip+1) = dedxmmy !(1) * ----------------------------------- CP_STORE(nplanes-ip+1,ibest)=idm if(DEBUG)print*,'%%%% included couple ',idm $ ,' (norm.dist.= ',distmin,', cut ',clinc,' )' goto 133 !next plane endif * ================================================ * STEP 2 >>>>>>> try to include a single cluster * either from a couple or single * ================================================ c if(DEBUG)print*,'>>>> try to include a new cluster' distmin=1000000. xmm_A = 0. !--------------------------- ymm_A = 0. ! init variables that zmm_A = 0. ! define the SINGLET xmm_B = 0. ! ymm_B = 0. ! zmm_B = 0. ! rxmm = 0. ! rymm = 0. ! dedxmmx = 0. !(1) dedxmmy = 0. !(1) iclm=0 !--------------------------- distance=0. *----- clusters inside couples ------------------------------------- do icp=1,ncp_plane(ip) !loop on cluster inside couples icx=clx(ip,icp) icy=cly(ip,icp) id=id_cp(ip,icp,ist) if(LADDER(icx).ne.nldt)goto 11882 !if the ladder number does not match * !jump to the next couple *----- try cluster x ----------------------------------------------- if(cl_used(icx).eq.1)goto 11881 !if the X cluster is already used * !jump to the Y cluster call xyz_PAM(icx,0,ist, c $ 'ETA2','ETA2', $ PFA,PFA, $ AXV_STORE(nplanes-ip+1,ibest),0.) distance = distance_to(XP,YP) c if(DEBUG)print*,'normalized distance ',distance if(DEBUG)print*,'( cl-X ',icx $ ,' in cp ',id,' ) normalized distance ',distance if(distance.lt.distmin)then xmm_A = xPAM_A ymm_A = yPAM_A zmm_A = zPAM_A xmm_B = xPAM_B ymm_B = yPAM_B zmm_B = zPAM_B rxmm = resxPAM rymm = resyPAM distmin = distance iclm = icx c dedxmm = dedx(icx) !(1) dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2) dedxmmy = 0. !(1) endif 11881 continue *----- try cluster y ----------------------------------------------- if(cl_used(icy).eq.1)goto 11882 !if the Y cluster is already used * !jump to the next couple call xyz_PAM(0,icy,ist, c $ 'ETA2','ETA2', $ PFA,PFA, $ 0.,AYV_STORE(nplanes-ip+1,ibest)) distance = distance_to(XP,YP) if(DEBUG)print*,'( cl-Y ',icy $ ,' in cp ',id,' ) normalized distance ',distance if(distance.lt.distmin)then xmm_A = xPAM_A ymm_A = yPAM_A zmm_A = zPAM_A xmm_B = xPAM_B ymm_B = yPAM_B zmm_B = zPAM_B rxmm = resxPAM rymm = resyPAM distmin = distance iclm = icy c dedxmm = dedx(icy) !(1) dedxmmx = 0. !(1) dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2) endif 11882 continue enddo !end loop on cluster inside couples *----- single clusters ----------------------------------------------- do ic=1,ncls(ip) !loop on single clusters icl=cls(ip,ic) if(cl_used(icl).eq.1.or. !if the cluster is already used $ LADDER(icl).ne.nldt.or. !or the ladder number does not match $ .false.)goto 18882 !jump to the next singlet if(mod(VIEW(icl),2).eq.0)then!<---- X view call xyz_PAM(icl,0,ist, c $ 'ETA2','ETA2', $ PFA,PFA, $ AXV_STORE(nplanes-ip+1,ibest),0.) else !<---- Y view call xyz_PAM(0,icl,ist, c $ 'ETA2','ETA2', $ PFA,PFA, $ 0.,AYV_STORE(nplanes-ip+1,ibest)) endif distance = distance_to(XP,YP) if(DEBUG)print*,'( cl-s ',icl $ ,' ) normalized distance ',distance if(distance.lt.distmin)then xmm_A = xPAM_A ymm_A = yPAM_A zmm_A = zPAM_A xmm_B = xPAM_B ymm_B = yPAM_B zmm_B = zPAM_B rxmm = resxPAM rymm = resyPAM distmin = distance iclm = icl c dedxmm = dedx(icl) !(1) if(mod(VIEW(icl),2).eq.0)then !<---- X view dedxmmx = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2) dedxmmy = 0. !(1) else !<---- Y view dedxmmx = 0. !(1) dedxmmy = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2) endif endif 18882 continue enddo !end loop on single clusters if(distmin.le.clinc)then CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<< * ---------------------------- if(mod(VIEW(iclm),2).eq.0)then XGOOD(nplanes-ip+1)=1. resx(nplanes-ip+1)=rxmm if(DEBUG)print*,'%%%% included X-cl ',iclm $ ,' ( norm.dist.= ',distmin,', cut ',clinc,' )' else YGOOD(nplanes-ip+1)=1. resy(nplanes-ip+1)=rymm if(DEBUG)print*,'%%%% included Y-cl ',iclm $ ,' ( norm.dist.= ',distmin,', cut ',clinc,' )' endif * ---------------------------- xm_A(nplanes-ip+1) = xmm_A ym_A(nplanes-ip+1) = ymm_A xm_B(nplanes-ip+1) = xmm_B ym_B(nplanes-ip+1) = ymm_B zm(nplanes-ip+1) = (zmm_A+zmm_B)/2. c dedxtrk(nplanes-ip+1) = dedxmm !<<< !(1) dedxtrk_x(nplanes-ip+1) = dedxmmx !<<< !(1) dedxtrk_y(nplanes-ip+1) = dedxmmy !<<< !(1) * ---------------------------- endif endif 133 continue enddo !end loop on planes return end *************************************************** * * * * * * * * * * * * ************************************************** subroutine clean_XYclouds(ibest,iflag) include '../common/commontracker.f' include '../common/common_momanhough.f' include '../common/momanhough_init.f' c include '../common/calib.f' c include '../common/level1.f' logical DEBUG common/dbg/DEBUG do ip=1,nplanes !loop on planes id=CP_STORE(nplanes-ip+1,ibest) icl=CLS_STORE(nplanes-ip+1,ibest) if(id.ne.0.or.icl.ne.0)then if(id.ne.0)then iclx=clx(ip,icp_cp(id)) icly=cly(ip,icp_cp(id)) cl_used(iclx)=1 !tag used clusters cl_used(icly)=1 !tag used clusters elseif(icl.ne.0)then cl_used(icl)=1 !tag used clusters endif c if(DEBUG)then c print*,ip,' <<< ',id c endif * ----------------------------- * remove the couple from clouds * remove also vitual couples containing the * selected clusters * ----------------------------- do icp=1,ncp_plane(ip) if( $ clx(ip,icp).eq.iclx $ .or. $ clx(ip,icp).eq.icl $ .or. $ cly(ip,icp).eq.icly $ .or. $ cly(ip,icp).eq.icl $ )then id=id_cp(ip,icp,1) if(DEBUG)then print*,ip,' <<< cp ',id $ ,' ( cl-x ' $ ,clx(ip,icp) $ ,' cl-y ' $ ,cly(ip,icp),' ) --> removed' endif * ----------------------------- * remove the couple from clouds do iyz=1,nclouds_yz if(cpcloud_yz(iyz,abs(id)).ne.0)then ptcloud_yz(iyz)=ptcloud_yz(iyz)-1 cpcloud_yz(iyz,abs(id))=0 endif enddo do ixz=1,nclouds_xz if(cpcloud_xz(ixz,abs(id)).ne.0)then ptcloud_xz(ixz)=ptcloud_xz(ixz)-1 cpcloud_xz(ixz,abs(id))=0 endif enddo * ----------------------------- endif enddo endif enddo !end loop on planes return end c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * *** c$$$ real function fbad_cog(ncog,ic) c$$$ c$$$ c$$$ include '../common/commontracker.f' c$$$ include '../common/level1.f' c$$$ include '../common/calib.f' c$$$ c$$$* --> signal of the central strip c$$$ sc = CLSIGNAL(INDMAX(ic)) !center c$$$ c$$$* signal of adjacent strips c$$$* --> left c$$$ sl1 = 0 !left 1 c$$$ if( c$$$ $ (INDMAX(ic)-1).ge.INDSTART(ic) c$$$ $ ) c$$$ $ sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1)) c$$$ c$$$ sl2 = 0 !left 2 c$$$ if( c$$$ $ (INDMAX(ic)-2).ge.INDSTART(ic) c$$$ $ ) c$$$ $ sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2)) c$$$ c$$$* --> right c$$$ sr1 = 0 !right 1 c$$$ if( c$$$ $ (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1)) c$$$ $ .or. c$$$ $ (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH) c$$$ $ ) c$$$ $ sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1)) c$$$ c$$$ sr2 = 0 !right 2 c$$$ if( c$$$ $ (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1)) c$$$ $ .or. c$$$ $ (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH) c$$$ $ ) c$$$ $ sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2)) c$$$ c$$$ c$$$ if(mod(int(VIEW(ic)),2).eq.1)then !Y-view c$$$ f = 4. c$$$ si = 8.4 c$$$ else !X-view c$$$ f = 6. c$$$ si = 3.9 c$$$ endif c$$$ c$$$ fbad_cog = 1. c$$$ f0 = 1 c$$$ f1 = 1 c$$$ f2 = 1 c$$$ f3 = 1 c$$$ if(sl1.gt.sr1.and.sl1.gt.0.)then c$$$ c$$$ if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)) ).eq.0)f0=f c$$$ if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f1=f c$$$c if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f3=f c$$$ c$$$ if(ncog.eq.2.and.sl1.ne.0)then c$$$ fbad_cog = (f1**2*sc**2/sl1**2+f0**2)/(sc**2/sl1**2+1.) c$$$ elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then c$$$ fbad_cog = 1. c$$$ elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then c$$$ fbad_cog = 1. c$$$ else c$$$ fbad_cog = 1. c$$$ endif c$$$ c$$$ elseif(sl1.le.sr1.and.sr1.gt.0.)then c$$$ c$$$ c$$$ if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)) ).eq.0)f0=f c$$$ if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f1=f c$$$c if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f3=f c$$$ c$$$ if(ncog.eq.2.and.sr1.ne.0)then c$$$ fbad_cog = (f1**2*sc**2/sr1**2+f0**2)/(sc**2/sr1**2+1.) c$$$ elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then c$$$ fbad_cog = 1. c$$$ elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then c$$$ fbad_cog = 1. c$$$ else c$$$ fbad_cog = 1. c$$$ endif c$$$ c$$$ endif c$$$ c$$$ fbad_cog = sqrt(fbad_cog) c$$$ c$$$ return c$$$ end c$$$