************************************************************ * The following subroutines * - track_finding >> hough transform * - track_fitting >> bob golden fitting * all the procedures to create LEVEL2 data, starting from LEVEL1 data. * * * * (This subroutine and all the dependent subroutines * will be included in the flight software) ************************************************************ subroutine track_finding(iflag) include 'commontracker.f' include 'level1.f' include 'common_momanhough.f' include 'common_mech.f' include 'common_xyzPAM.f' include 'common_mini_2.f' include 'calib.f' include 'level2.f' c print*,'======================================================' c$$$ do ic=1,NCLSTR1 c$$$ if(.false. c$$$ $ .or.nsatstrips(ic).gt.0 c$$$c $ .or.nbadstrips(0,ic).gt.0 c$$$c $ .or.nbadstrips(4,ic).gt.0 c$$$c $ .or.nbadstrips(3,ic).gt.0 c$$$ $ .or..false.)then c$$$ print*,'--- cl-',ic,' ------------------------' c$$$ istart = INDSTART(IC) c$$$ istop = TOTCLLENGTH c$$$ if(ic.lt.NCLSTR1)istop=INDSTART(IC+1)-1 c$$$ print*,'ADC ',(CLADC(i),i=istart,istop) c$$$ print*,'s/n ',(CLSIGNAL(i)/CLSIGMA(i),i=istart,istop) c$$$ print*,'sgnl ',(CLSIGNAL(i),i=istart,istop) c$$$ print*,'strip ',(i-INDMAX(ic),i=istart,istop) c$$$ print*,'view ',VIEW(ic) c$$$ print*,'maxs ',MAXS(ic) c$$$ print*,'COG4 ',cog(4,ic) c$$$ ff = fbad_cog(4,ic) c$$$ print*,'fbad ',ff c$$$ print*,(CLBAD(i),i=istart,istop) c$$$ bb=nbadstrips(0,ic) c$$$ print*,'#BAD (tot)',bb c$$$ bb=nbadstrips(4,ic) c$$$ print*,'#BAD (4)',bb c$$$ bb=nbadstrips(3,ic) c$$$ print*,'#BAD (3)',bb c$$$ ss=nsatstrips(ic) c$$$ print*,'#saturated ',ss c$$$ endif c$$$ enddo *------------------------------------------------------------------------------- * STEP 1 *------------------------------------------------------------------------------- * X-Y cluster association * * Clusters are associated to form COUPLES * Clusters not associated in any couple are called SINGLETS * * Track identification (Hough transform) and fitting is first done on couples. * Hence singlets are possibly added to the track. * * Variables assigned by the routine "cl_to_couples" are those in the * common blocks: * - common/clusters/cl_good * - common/couples/clx,cly,ncp_plane,ncp_tot,cp_useds1,cp_useds2 * - common/singlets/ncls,cls,cl_single *------------------------------------------------------------------------------- *------------------------------------------------------------------------------- call cl_to_couples(iflag) if(iflag.eq.1)then !bad event goto 880 !go to next event endif if(ncp_tot.eq.0)goto 880 !go to next event *----------------------------------------------------- *----------------------------------------------------- * HOUGH TRASFORM *----------------------------------------------------- *----------------------------------------------------- *------------------------------------------------------------------------------- * STEP 2 *------------------------------------------------------------------------------- * * Association of couples to form * - DOUBLETS in YZ view * - TRIPLETS in XZ view * * Variables assigned by the routine "cp_to_doubtrip" are those in the * common blocks: * - common/hough_param/ * $ alfayz1, !Y0 * $ alfayz2, !tg theta-yz * $ alfaxz1, !X0 * $ alfaxz2, !tg theta-xz * $ alfaxz3 !1/r * - common/doublets/ndblt,cpyz1,cpyz2 * - common/triplets/ntrpt,cpxz1,cpxz2,cpxz3 *------------------------------------------------------------------------------- *------------------------------------------------------------------------------- call cp_to_doubtrip(iflag) if(iflag.eq.1)then !bad event goto 880 !go to next event endif if(ndblt.eq.0.or.ntrpt.eq.0)goto 880 !go to next event *------------------------------------------------------------------------------- * STEP 3 *------------------------------------------------------------------------------- * * Classification of doublets and triplets to form CLOUDS, * according to distance in parameter space. * * cloud = cluster of points (doublets/triplets) in parameter space * * * * Variables assigned by the routine "doub_to_YZcloud" are those in the * common blocks: * - common/clouds_yz/ * $ nclouds_yz * $ ,alfayz1_av,alfayz2_av * $ ,ptcloud_yz,db_cloud,cpcloud_yz * * Variables assigned by the routine "trip_to_XZcloud" are those in the * common blocks: * common/clouds_xz/ * $ nclouds_xz xz2_av,alfaxz3_av * $ ,ptcloud_xz,tr_cloud,cpcloud_xz *------------------------------------------------------------------------------- *------------------------------------------------------------------------------- * count number of hit planes planehit=0 do np=1,nplanes if(ncp_plane(np).ne.0)then planehit=planehit+1 endif enddo if(planehit.lt.3) goto 880 ! exit nptxz_min=x_min_start nplxz_min=x_min_start nptyz_min=y_min_start nplyz_min=y_min_start cutdistyz=cutystart cutdistxz=cutxstart 878 continue call doub_to_YZcloud(iflag) if(iflag.eq.1)then !bad event goto 880 !fill ntp and go to next event endif * ------------------------------------------------ * first try to release the tolerance * ------------------------------------------------ if(nclouds_yz.eq.0.and.cutdistyz.lt.maxcuty)then if(cutdistyz.lt.maxcuty/2)then cutdistyz=cutdistyz+cutystep else cutdistyz=cutdistyz+(3*cutystep) endif if(DEBUG.EQ.1)print*,'try again - cutdistyz ',cutdistyz goto 878 endif * ------------------------------------------------ * hence reduce the minimum number of plane * ------------------------------------------------ if(nclouds_yz.eq.0.and.nplyz_min.gt.3)then nplyz_min=nplyz_min-1 if(DEBUG.EQ.1)print*,'try again - nplyz_min ',nplyz_min goto 878 endif if(nclouds_yz.eq.0)goto 880 !go to next event ccc if(planehit.eq.3) goto 881 879 continue call trip_to_XZcloud(iflag) if(iflag.eq.1)then !bad event goto 880 !fill ntp and go to next event endif * ------------------------------------------------ * first try to release the tolerance * ------------------------------------------------ if(nclouds_xz.eq.0.and.cutdistxz.lt.maxcutx)then cutdistxz=cutdistxz+cutxstep if(DEBUG.EQ.1)print*,'try again - cutdistxz ',cutdistxz goto 879 endif * ------------------------------------------------ * hence reduce the minimum number of plane * ------------------------------------------------ if(nclouds_xz.eq.0.and.nplxz_min.gt.3)then nplxz_min=nplxz_min-1 if(DEBUG.EQ.1)print*,'try again - nplxz_min ',nplxz_min goto 879 endif if(nclouds_xz.eq.0)goto 880 !go to next event c$$$ 881 continue c$$$* if there is at least three planes on the Y view decreases cuts on X view c$$$ if(nclouds_xz.eq.0.and.nclouds_yz.gt.0.and. c$$$ $ nplxz_min.ne.y_min_start)then c$$$ nptxz_min=x_min_step c$$$ nplxz_min=x_min_start-x_min_step c$$$ goto 879 c$$$ endif 880 return end ************************************************************ subroutine track_fitting(iflag) include 'commontracker.f' include 'level1.f' include 'common_momanhough.f' include 'common_mech.f' include 'common_xyzPAM.f' include 'common_mini_2.f' include 'calib.f' include 'level2.f' c include 'momanhough_init.f' logical FIMAGE ! real trackimage(NTRACKSMAX) real*8 AL_GUESS(5) *------------------------------------------------------------------------------- * STEP 4 (ITERATED until any other physical track isn't found) *------------------------------------------------------------------------------- * * YZ and XZ clouds are combined in order to obtain the initial guess * of the candidate-track parameters. * A minimum number of matching couples between YZ and XZ clouds is required. * * A TRACK CANDIDATE is defined by * - the couples resulting from the INTERSECTION of the two clouds, and * - the associated track parameters (evaluated by performing a zero-order * track fitting) * * The NTRACKS candidate-track parameters are stored in common block: * * - common/track_candidates/NTRACKS,AL_STORE * $ ,XV_STORE,YV_STORE,ZV_STORE * $ ,XM_STORE,YM_STORE,ZM_STORE * $ ,RESX_STORE,RESY_STORE * $ ,AXV_STORE,AYV_STORE * $ ,XGOOD_STORE,YGOOD_STORE * $ ,CP_STORE,RCHI2_STORE * *------------------------------------------------------------------------------- *------------------------------------------------------------------------------- ccc ntrk=0 !counter of identified physical tracks c11111 continue !<<<<<<< come here when performing a new search continue !<<<<<<< come here when performing a new search if(nclouds_xz.eq.0)goto 880 !go to next event if(nclouds_yz.eq.0)goto 880 !go to next event c iflag=0 call clouds_to_ctrack(iflag) if(iflag.eq.1)then !no candidate tracks found goto 880 !fill ntp and go to next event endif if(ntracks.eq.0)goto 880 !go to next event FIMAGE=.false. !processing best track (not track image) ibest=0 !best track among candidates iimage=0 !track image * ------------- select the best track ------------- c$$$ rchi2best=1000000000. c$$$ do i=1,ntracks c$$$ if(RCHI2_STORE(i).lt.rchi2best.and. c$$$ $ RCHI2_STORE(i).gt.0)then c$$$ ibest=i c$$$ rchi2best=RCHI2_STORE(i) c$$$ endif c$$$ enddo c$$$ if(ibest.eq.0)goto 880 !>> no good candidates * ------------------------------------------------------- * order track-candidates according to: * 1st) decreasing n.points * 2nd) increasing chi**2 * ------------------------------------------------------- rchi2best=1000000000. ndofbest=0 do i=1,ntracks ndof=0 do ii=1,nplanes ndof=ndof $ +int(xgood_store(ii,i)) $ +int(ygood_store(ii,i)) enddo if(ndof.gt.ndofbest)then ibest=i rchi2best=RCHI2_STORE(i) ndofbest=ndof elseif(ndof.eq.ndofbest)then if(RCHI2_STORE(i).lt.rchi2best.and. $ RCHI2_STORE(i).gt.0)then ibest=i rchi2best=RCHI2_STORE(i) ndofbest=ndof endif endif enddo if(ibest.eq.0)goto 880 !>> no good candidates *------------------------------------------------------------------------------- * The best track candidate (ibest) is selected and a new fitting is performed. * Previous to this, the track is refined by: * - possibly adding new COUPLES or SINGLETS from the missing planes * - evaluating the coordinates with improved PFAs * ( angle-dependent ETA algorithms ) *------------------------------------------------------------------------------- 1212 continue !<<<<< come here to fit track-image if(.not.FIMAGE)then !processing best candidate icand=ibest else !processing image icand=iimage iimage=0 endif if(icand.eq.0)then if(VERBOSE.EQ.1)then print*,'HAI FATTO UN CASINO!!!!!! icand = ',icand $ ,ibest,iimage endif return endif * *-*-*-*-*-*-*-*-*-*-*-*-*-*-* call refine_track(icand) * *-*-*-*-*-*-*-*-*-*-*-*-*-*-* * ********************************************************** * ************************** FIT *** FIT *** FIT *** FIT *** * ********************************************************** call guess() do i=1,5 AL_GUESS(i)=AL(i) enddo do i=1,5 AL(i)=dble(AL_STORE(i,icand)) enddo IDCAND = icand !fitted track-candidate ifail=0 !error flag in chi2 computation jstep=0 !# minimization steps iprint=0 c if(DEBUG.EQ.1)iprint=1 if(VERBOSE.EQ.1)iprint=1 if(DEBUG.EQ.1)iprint=2 call mini2(jstep,ifail,iprint) cc if(ifail.ne.0) then if(ifail.ne.0.or.CHI2.ne.CHI2) then !new if(CHI2.ne.CHI2)CHI2=-9999. !new if(VERBOSE.EQ.1)then print *, $ '*** MINIMIZATION FAILURE *** (after refinement) ' $ ,iev endif endif if(DEBUG.EQ.1)then print*,'----------------------------- improved track coord' 22222 format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5) do ip=1,6 write(*,22222)ip,zm(ip),xm(ip),ym(ip) $ ,xm_A(ip),ym_A(ip),xm_B(ip),ym_B(ip) $ ,xgood(ip),ygood(ip),resx(ip),resy(ip) enddo endif c rchi2=chi2/dble(ndof) if(DEBUG.EQ.1)then print*,' ' print*,'****** SELECTED TRACK *************' print*,'# R. chi2 RIG' print*,' --- ',chi2,' --- ' $ ,1./abs(AL(5)) print*,'***********************************' endif * ********************************************************** * ************************** FIT *** FIT *** FIT *** FIT *** * ********************************************************** * ------------- search if the track has an IMAGE ------------- * ------------- (also this is stored ) ------------- if(FIMAGE)goto 122 !>>> jump! (this is already an image) * ----------------------------------------------------- * first check if the track is ambiguous * ----------------------------------------------------- * (modified on august 2007 by ElenaV) is1=0 do ip=1,NPLANES if(SENSOR_STORE(ip,icand).ne.0)then is1=SENSOR_STORE(ip,icand) if(ip.eq.6)is1=3-is1 !last plane inverted endif enddo if(is1.eq.0)then if(WARNING.EQ.1)print*,'** WARNING ** sensor=0' goto 122 !jump endif do ip=1,NPLANES is2 = SENSOR_STORE(ip,icand) !sensor if(ip.eq.6.and.is2.ne.0)is2=3-is2 !last plane inverted if( $ (is1.ne.is2.and.is2.ne.0) $ )goto 122 !jump (this track cannot have an image) enddo if(DEBUG.eq.1)print*,' >>> ambiguous track! ' * --------------------------------------------------------------- * take the candidate with the greatest number of matching couples * if more than one satisfies the requirement, * choose the one with more points and lower chi2 * --------------------------------------------------------------- * count the number of matching couples ncpmax = 0 iimage = 0 !id of candidate with better matching do i=1,ntracks ncp=0 do ip=1,nplanes if(CP_STORE(nplanes-ip+1,icand).ne.0)then if( $ CP_STORE(nplanes-ip+1,i).ne.0 $ .and. $ CP_STORE(nplanes-ip+1,icand).eq. $ -1*CP_STORE(nplanes-ip+1,i) $ )then ncp=ncp+1 !ok elseif( $ CP_STORE(nplanes-ip+1,i).ne.0 $ .and. $ CP_STORE(nplanes-ip+1,icand).ne. $ -1*CP_STORE(nplanes-ip+1,i) $ )then ncp = 0 goto 117 !it is not an image candidate else endif endif enddo 117 continue trackimage(i)=ncp !number of matching couples if(ncp>ncpmax)then ncpmax=ncp iimage=i endif enddo * check if there are more than one image candidates ngood=0 do i=1,ntracks if( ncpmax.ne.0.and.trackimage(i).eq.ncpmax )ngood=ngood+1 enddo if(DEBUG.eq.1)print*,' n.image-candidates : ',ngood * if there are, choose the best one if(ngood.gt.0)then * ------------------------------------------------------- * order track-candidates according to: * 1st) decreasing n.points * 2nd) increasing chi**2 * ------------------------------------------------------- rchi2best=1000000000. ndofbest=0 do i=1,ntracks if( trackimage(i).eq.ncpmax )then ndof=0 do ii=1,nplanes ndof=ndof $ +int(xgood_store(ii,i)) $ +int(ygood_store(ii,i)) enddo if(ndof.gt.ndofbest)then iimage=i rchi2best=RCHI2_STORE(i) ndofbest=ndof elseif(ndof.eq.ndofbest)then if(RCHI2_STORE(i).lt.rchi2best.and. $ RCHI2_STORE(i).gt.0)then iimage=i rchi2best=RCHI2_STORE(i) ndofbest=ndof endif endif endif enddo if(DEBUG.EQ.1)then print*,'Track candidate ',iimage $ ,' >>> TRACK IMAGE >>> of' $ ,ibest endif endif 122 continue * --- and store the results -------------------------------- if(ntrk.eq.NTRKMAX)then if(verbose.eq.1) $ print*, $ '** warning ** number of identified '// $ 'tracks exceeds vector dimension ' $ ,'( ',NTRKMAX,' )' cc good2=.false. goto 880 !fill ntp and go to next event endif ntrk = ntrk + 1 !counter of found tracks if(.not.FIMAGE $ .and.iimage.eq.0) image(ntrk)= 0 if(.not.FIMAGE $ .and.iimage.ne.0)image(ntrk)=ntrk+1 !this is the image of the next if(FIMAGE) image(ntrk)=ntrk-1 !this is the image of the previous call fill_level2_tracks(ntrk) !==> good2=.true. c$$$ if(ntrk.eq.NTRKMAX)then c$$$ if(verbose.eq.1) c$$$ $ print*, c$$$ $ '** warning ** number of identified '// c$$$ $ 'tracks exceeds vector dimension ' c$$$ $ ,'( ',NTRKMAX,' )' c$$$cc good2=.false. c$$$ goto 880 !fill ntp and go to next event c$$$ endif if(iimage.ne.0)then FIMAGE=.true. ! goto 1212 !>>> fit image-track endif 880 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 * bfx - x-component of magnetci field * bfy - y-component of magnetic field * * --------- 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_xyzPAM.f) * * subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy) include 'commontracker.f' include 'level1.f' include 'calib.f' include 'common_align.f' include 'common_mech.f' include 'common_xyzPAM.f' integer icx,icy !X-Y cluster ID integer sensor integer viewx,viewy character*4 PFAx,PFAy !PFA to be used real ax,ay !X-Y geometric angle real angx,angy !X-Y effective angle real bfx,bfy !X-Y b-field components real stripx,stripy double precision xi,yi,zi double precision xi_A,yi_A,zi_A double precision xi_B,yi_B,zi_B double precision xrt,yrt,zrt double precision xrt_A,yrt_A,zrt_A double precision xrt_B,yrt_B,zrt_B parameter (ndivx=30) resxPAM = 0 resyPAM = 0 xPAM = 0.D0 yPAM = 0.D0 zPAM = 0.D0 xPAM_A = 0.D0 yPAM_A = 0.D0 zPAM_A = 0.D0 xPAM_B = 0.D0 yPAM_B = 0.D0 zPAM_B = 0.D0 if(sensor.lt.1.or.sensor.gt.2)then print*,'xyz_PAM ***ERROR*** wrong input ' print*,'sensor ',sensor icx=0 icy=0 endif * ----------------- * CLUSTER X * ----------------- if(icx.ne.0)then viewx = VIEW(icx) nldx = nld(MAXS(icx),VIEW(icx)) nplx = npl(VIEW(icx)) c resxPAM = RESXAV stripx = float(MAXS(icx)) if( $ viewx.lt.1.or. $ viewx.gt.12.or. $ nldx.lt.1.or. $ nldx.gt.3.or. $ stripx.lt.1.or. $ stripx.gt.3072.or. $ .false.)then print*,'xyz_PAM ***ERROR*** wrong input ' print*,'icx ',icx,'view ',viewx,'nld ',nldx,'strip ',stripx icx = 0 goto 10 endif * -------------------------- * magnetic-field corrections * -------------------------- stripx = stripx + fieldcorr(viewx,bfy) angx = effectiveangle(ax,viewx,bfy) call applypfa(PFAx,icx,angx,corr,res) stripx = stripx + corr resxPAM = res 10 continue endif * ----------------- * CLUSTER Y * ----------------- if(icy.ne.0)then viewy = VIEW(icy) nldy = nld(MAXS(icy),VIEW(icy)) nply = npl(VIEW(icy)) c resyPAM = RESYAV stripy = float(MAXS(icy)) if( $ viewy.lt.1.or. $ viewy.gt.12.or. $ nldy.lt.1.or. $ nldy.gt.3.or. $ stripy.lt.1.or. $ stripy.gt.3072.or. $ .false.)then print*,'xyz_PAM ***ERROR*** wrong input ' print*,'icy ',icy,'view ',viewy,'nld ',nldy,'strip ',stripy icy = 0 goto 20 endif if(icx.ne.0.and.(nply.ne.nplx.or.nldy.ne.nldx))then if(DEBUG.EQ.1) then print*,'xyz_PAM ***ERROR*** invalid cluster couple!!! ' $ ,icx,icy endif goto 100 endif * -------------------------- * magnetic-field corrections * -------------------------- stripy = stripy + fieldcorr(viewy,bfx) angy = effectiveangle(ay,viewy,bfx) call applypfa(PFAy,icy,angy,corr,res) stripy = stripy + corr resyPAM = res 20 continue 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------------------------------------------------------------------------ if(((mod(int(stripx+0.5)-1,1024)+1).le.3) $ .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips... if(DEBUG.EQ.1) then print*,'xyz_PAM (couple):', $ ' WARNING: false X strip: strip ',stripx endif endif xi = dcoordsi(stripx,viewx) yi = dcoordsi(stripy,viewy) zi = 0.D0 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.D0 yPAM_A = 0.D0 zPAM_A = 0.D0 xPAM_B = 0.D0 yPAM_B = 0.D0 zPAM_B = 0.D0 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 xi = dcoordsi(0.5*(nstrips+1),viewx) !sensor center yi = dcoordsi(stripy,viewy) zi = 0.D0 xi_A = edgeY_d - SiDimX/2 yi_A = yi zi_A = 0. xi_B = SiDimX/2 - edgeY_u yi_B = yi zi_B = 0. elseif(icx.ne.0)then c=========================================================== C X-SINGLET C=========================================================== nply = nplx nldy = nldx viewy = viewx - 1 if(((mod(int(stripx+0.5)-1,1024)+1).le.3) $ .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips... if(DEBUG.EQ.1) then print*,'xyz_PAM (X-singlet):', $ ' WARNING: false X strip: strip ',stripx endif endif xi = dcoordsi(stripx,viewx) yi = dcoordsi(0.5*(nstrips+1),viewy) !sensor center zi = 0.D0 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 else if(DEBUG.EQ.1) then print *,'routine xyz_PAM ---> not properly used !!!' print *,'icx = ',icx print *,'icy = ',icy endif 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) 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 c$$$ xPAM = 0.D0 c$$$ yPAM = 0.D0 c$$$ zPAM = 0.D0 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 else if(DEBUG.EQ.1) then print *,'routine xyz_PAM ---> not properly used !!!' print *,'icx = ',icx print *,'icy = ',icy endif endif 100 continue end ************************************************************************ * Call xyz_PAM subroutine with default PFA and fill the mini2 common. * (done to be called from c/c++) ************************************************************************ subroutine xyzpam(ip,icx,icy,lad,sensor,ax,ay,bfx,bfy) include 'commontracker.f' include 'level1.f' include 'common_mini_2.f' include 'common_xyzPAM.f' include 'common_mech.f' include 'calib.f' * flag to chose PFA c$$$ character*10 PFA c$$$ common/FINALPFA/PFA integer icx,icy !X-Y cluster ID integer sensor character*4 PFAx,PFAy !PFA to be used real ax,ay !X-Y geometric angle real bfx,bfy !X-Y b-field components ipx=0 ipy=0 c$$$ PFAx = 'COG4'!PFA c$$$ PFAy = 'COG4'!PFA if(icx.gt.nclstr1.or.icy.gt.nclstr1)then print*,'xyzpam: ***WARNING*** clusters ',icx,icy $ ,' do not exists (n.clusters=',nclstr1,')' icx = -1*icx icy = -1*icy return endif call idtoc(pfaid,PFAx) call idtoc(pfaid,PFAy) if(icx.ne.0.and.icy.ne.0)then ipx=npl(VIEW(icx)) ipy=npl(VIEW(icy)) if( (nplanes-ipx+1).ne.ip )then print*,'xyzpam: ***WARNING*** cluster ',icx $ ,' does not belong to plane: ',ip icx = -1*icx return endif if( (nplanes-ipy+1).ne.ip )then print*,'xyzpam: ***WARNING*** cluster ',icy $ ,' does not belong to plane: ',ip icy = -1*icy return endif call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy) xgood(ip) = 1. ygood(ip) = 1. resx(ip) = resxPAM resy(ip) = resyPAM xm(ip) = xPAM ym(ip) = yPAM zm(ip) = zPAM xm_A(ip) = 0.D0 ym_A(ip) = 0.D0 xm_B(ip) = 0.D0 ym_B(ip) = 0.D0 c zv(ip) = zPAM elseif(icx.eq.0.and.icy.ne.0)then ipy=npl(VIEW(icy)) if( (nplanes-ipy+1).ne.ip )then print*,'xyzpam: ***WARNING*** cluster ',icy $ ,' does not belong to plane: ',ip icy = -1*icy return endif call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy) xgood(ip) = 0. ygood(ip) = 1. resx(ip) = 1000. resy(ip) = resyPAM c$$$ xm(ip) = -100. c$$$ ym(ip) = -100. c$$$ zm(ip) = (zPAM_A+zPAM_B)/2. xm(ip) = xPAM ym(ip) = yPAM zm(ip) = zPAM xm_A(ip) = xPAM_A ym_A(ip) = yPAM_A xm_B(ip) = xPAM_B ym_B(ip) = yPAM_B c zv(ip) = (zPAM_A+zPAM_B)/2. elseif(icx.ne.0.and.icy.eq.0)then ipx=npl(VIEW(icx)) if( (nplanes-ipx+1).ne.ip )then print*,'xyzpam: ***WARNING*** cluster ',icx $ ,' does not belong to plane: ',ip icx = -1*icx return endif call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy) xgood(ip) = 1. ygood(ip) = 0. resx(ip) = resxPAM resy(ip) = 1000. c$$$ xm(ip) = -100. c$$$ ym(ip) = -100. c$$$ zm(ip) = (zPAM_A+zPAM_B)/2. xm(ip) = xPAM ym(ip) = yPAM zm(ip) = zPAM xm_A(ip) = xPAM_A ym_A(ip) = yPAM_A xm_B(ip) = xPAM_B ym_B(ip) = yPAM_B c zv(ip) = (zPAM_A+zPAM_B)/2. else il = 2 if(lad.ne.0)il=lad is = 1 if(sensor.ne.0)is=sensor xgood(ip) = 0. ygood(ip) = 0. resx(ip) = 1000. resy(ip) = 1000. xm(ip) = -100. ym(ip) = -100. zm(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4 xm_A(ip) = 0. ym_A(ip) = 0. xm_B(ip) = 0. ym_B(ip) = 0. c zv(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4 endif if(DEBUG.EQ.1)then 22222 format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5) write(*,22222)ip,zm(ip),xm(ip),ym(ip) $ ,xm_A(ip),ym_A(ip),xm_B(ip),ym_B(ip) $ ,xgood(ip),ygood(ip),resx(ip),resy(ip) endif 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(rXPP,rYPP) include 'common_xyzPAM.f' * ----------------------------------- * it computes the normalized distance * ( i.e. distance/resolution ) * ----------------------------------- real rXPP,rYPP double precision XPP,YPP double precision distance,RE double precision BETA,ALFA,xmi,ymi XPP=DBLE(rXPP) YPP=DBLE(rYPP) * ---------------------- if ( c$$$ + xPAM.eq.0.and. c$$$ + yPAM.eq.0.and. c$$$ + 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)!QUIQUI cc $ ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2 distance=dsqrt(distance) * ---------------------- elseif( c$$$ + xPAM.ne.0.and. c$$$ + yPAM.ne.0.and. c$$$ + 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))**2 !QUIQUI $ + $ ((yPAM-YPP))**2 c$$$ $ ((xPAM-XPP)/resxPAM)**2 c$$$ $ + c$$$ $ ((yPAM-YPP)/resyPAM)**2 distance=dsqrt(distance) else 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 'commontracker.f' include '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./ double precision yvvv,xvvv double precision xi,yi,zi double precision xrt,yrt,zrt real AA,BB double precision yvv(4),xvv(4) * tollerance to consider the track inside the sensitive area real ptoll data ptoll/150./ !um external nviewx,nviewy,dcoordsi,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 = dcoordsi(stripx,viewx) yi = dcoordsi(stripy,viewy) zi = 0.D0 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) 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 = REAL((yvv(iv1)-yvv(iv2))/(xvv(iv1)-xvv(iv2))) !EM GCC4.7 BB = REAL(yvv(iv1) - AA*xvv(iv1)) !EM GCC4.7 * 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 = REAL((xvv(iv1)-xvv(iv2))/(yvv(iv1)-yvv(iv2))) !EM GCC4.7 BB = REAL(xvv(iv1) - AA*yvv(iv1)) !EM GCC4.7 * 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 'commontracker.f' include 'level1.f' c include 'common_analysis.f' include '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 return end integer function icp_cp(id) * * given the couple id, * it returns the id number ON THE PLANE * include 'commontracker.f' include 'level1.f' c include 'common_analysis.f' include '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 'commontracker.f' include 'level1.f' c include 'calib.f' c include 'level1.f' c include 'common_analysis.f' include '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 cl_to_couples(iflag) include 'commontracker.f' include 'level1.f' include 'common_momanhough.f' c include 'momanhough_init.f' include 'calib.f' c include 'level1.f' * output flag * -------------- * 0 = good event * 1 = bad event * -------------- integer iflag integer badseed,badclx,badcly iflag = iflag if(DEBUG.EQ.1)print*,'cl_to_couples:' cc if(RECOVER_SINGLETS.and..not.SECOND_SEARCH)goto 80 * init variables 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 !all are single cl_good(icl) = 0 !all are bad enddo do iv=1,nviews ncl_view(iv) = 0 mask_view(iv) = 0 !all included enddo * count number of cluster per view do icl=1,nclstr1 ncl_view(VIEW(icl)) = ncl_view(VIEW(icl)) + 1 enddo * mask views with too many clusters do iv=1,nviews if( ncl_view(iv).gt. nclusterlimit)then c mask_view(iv) = 1 mask_view(iv) = mask_view(iv) + 2**0 if(DEBUG.EQ.1) $ print*,' * WARNING * cl_to_couple: n.clusters > ' $ ,nclusterlimit,' on view ', iv,' --> masked!' endif enddo * start association ncouples=0 do icx=1,nclstr1 !loop on cluster (X) if(mod(VIEW(icx),2).eq.1)goto 10 if(cl_used(icx).ne.0)goto 10 * ---------------------------------------------------- * jump masked views (X VIEW) * ---------------------------------------------------- if( mask_view(VIEW(icx)).ne.0 ) goto 10 * ---------------------------------------------------- * cut on charge (X VIEW) * ---------------------------------------------------- if(sgnl(icx)/mip(VIEW(icx),LADDER(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 badclx=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 badclx=badclx*ibad enddo * ---------------------------------------------------- * >>> eliminato il taglio sulle BAD <<< * ---------------------------------------------------- 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 if(cl_used(icx).ne.0)goto 20 * ---------------------------------------------------- * jump masked views (Y VIEW) * ---------------------------------------------------- if( mask_view(VIEW(icy)).ne.0 ) goto 20 * ---------------------------------------------------- * cut on charge (Y VIEW) * ---------------------------------------------------- if(sgnl(icy)/mip(VIEW(icy),LADDER(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 badcly=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)) badcly=badcly*ibad enddo * ---------------------------------------------------- * >>> eliminato il taglio sulle BAD <<< * ---------------------------------------------------- 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 * (modified to be applied only below saturation... obviously) if( .not.(sgnl(icy).gt.chsaty.and.sgnl(icx).gt.chsatx) $ .and. $ .not.(sgnl(icy).lt.chmipy.and.sgnl(icx).lt.chmipx) $ .and. $ (badclx.eq.1.and.badcly.eq.1) $ .and. $ .true.)then ddd=(sgnl(icy) $ -kch(nplx,nldx)*sgnl(icx)-cch(nplx,nldx)) ddd=ddd/sqrt(kch(nplx,nldx)**2+1) c cut = chcut * sch(nplx,nldx) sss=(kch(nplx,nldx)*sgnl(icy)+sgnl(icx) $ -kch(nplx,nldx)*cch(nplx,nldx)) sss=sss/sqrt(kch(nplx,nldx)**2+1) cut = chcut * (16 + sss/50.) if(abs(ddd).gt.cut)then goto 20 !charge not consistent endif endif if(ncp_plane(nplx).gt.ncouplemax)then if(verbose.eq.1)print*, $ '** warning ** number of identified '// $ 'couples on plane ',nplx, $ 'exceeds vector dimention ' $ ,'( ',ncouplemax,' ) --> masked!' c mask_view(nviewx(nplx)) = 2 c mask_view(nviewy(nply)) = 2 mask_view(nviewx(nplx))= mask_view(nviewx(nplx))+ 2**1 mask_view(nviewy(nply))= mask_view(nviewy(nply))+ 2**1 goto 10 endif * ------------------> COUPLE <------------------ 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 c 80 continue continue if(DEBUG.EQ.1)then print*,'clusters ',nclstr1 print*,'good ',(cl_good(i),i=1,nclstr1) print*,'used ',(cl_used(i),i=1,nclstr1) print*,'singlets',(cl_single(i),i=1,nclstr1) print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes) endif if(.not.RECOVER_SINGLETS)goto 81 * //////////////////////////////////////////////// * PATCH to recover events with less than 3 couples * //////////////////////////////////////////////// * loop over singlet and create "fake" couples * (with clx=0 or cly=0) * if(DEBUG.EQ.1) $ print*,'>>> Recover singlets ' $ ,'(creates fake couples) <<<' do icl=1,nclstr1 if( $ cl_single(icl).eq.1.and. $ cl_used(icl).eq.0.and. $ .true.)then * ---------------------------------------------------- * jump masked views * ---------------------------------------------------- if( mask_view(VIEW(icl)).ne.0 ) goto 21 if(mod(VIEW(icl),2).eq.0)then !=== X views * ---------------------------------------------------- * cut on charge (X VIEW) * ---------------------------------------------------- if(sgnl(icl).lt.dedx_x_min) goto 21 nplx=npl(VIEW(icl)) * ------------------> (FAKE) COUPLE <----------------- ncp_plane(nplx) = ncp_plane(nplx) + 1 clx(nplx,ncp_plane(nplx))=icl cly(nplx,ncp_plane(nplx))=0 c$$$ cl_single(icl)=0! I leave the cluster tagged as singlet!!! * ---------------------------------------------------- else !=== Y views * ---------------------------------------------------- * cut on charge (Y VIEW) * ---------------------------------------------------- if(sgnl(icl).lt.dedx_y_min) goto 21 nply=npl(VIEW(icl)) * ------------------> (FAKE) COUPLE <----------------- ncp_plane(nply) = ncp_plane(nply) + 1 clx(nply,ncp_plane(nply))=0 cly(nply,ncp_plane(nply))=icl c$$$ cl_single(icl)=0! I leave the cluster tagged as singlet!!! * ---------------------------------------------------- endif endif !end "single" condition 21 continue enddo !end loop over clusters if(DEBUG.EQ.1) $ print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes) 81 continue ncp_tot=0 do ip=1,NPLANES ncp_tot = ncp_tot + ncp_plane(ip) enddo if(DEBUG.EQ.1) $ print*,'n.couple tot: ',ncp_tot return end *************************************************** * * * * * * * * * * * * ************************************************** subroutine cp_to_doubtrip(iflag) include 'commontracker.f' include 'level1.f' include 'common_momanhough.f' include 'common_xyzPAM.f' include 'common_mini_2.f' include 'calib.f' * 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 * ----------------------------- if(DEBUG.EQ.1)print*,'cp_to_doubtrip:' * -------------------------------------------- * put a limit to the maximum number of couples * per plane, in order to apply hough transform * (couples recovered during track refinement) * -------------------------------------------- do ip=1,nplanes if(ncp_plane(ip).gt.ncouplelimit)then mask_view(nviewx(ip)) = mask_view(nviewx(ip)) + 2**7 mask_view(nviewy(ip)) = mask_view(nviewy(ip)) + 2**7 endif enddo ndblt=0 !number of doublets ntrpt=0 !number of triplets do ip1=1,(nplanes-1) !loop on planes - COPPIA 1 c$$$ print*,'(1) ip ',ip1 c$$$ $ ,mask_view(nviewx(ip1)) c$$$ $ ,mask_view(nviewy(ip1)) if( mask_view(nviewx(ip1)).ne.0 .or. $ mask_view(nviewy(ip1)).ne.0 )goto 10 !skip plane 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$$$ print*,'(1) ip ',ip1,' icp ',icp1 c call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1) c call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1) call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.,0.,0.) xm1=REAL(xPAM) !EM GCC4.7 ym1=REAL(yPAM) !EM GCC4.7 zm1=REAL(zPAM) !EM GCC4.7 do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2 c$$$ print*,'(2) ip ',ip2 c$$$ $ ,mask_view(nviewx(ip2)) c$$$ $ ,mask_view(nviewy(ip2)) if( mask_view(nviewx(ip2)).ne.0 .or. $ mask_view(nviewy(ip2)).ne.0 )goto 20 !skip plane 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$$$ print*,'(2) ip ',ip2,' icp ',icp2 c call xyz_PAM c $ (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1) c call xyz_PAM c $ (icx2,icy2,is2,PFAdef,PFAdef,0.,0.) !(1) call xyz_PAM $ (icx2,icy2,is2,PFAdef,PFAdef,0.,0.,0.,0.) xm2=REAL(xPAM) !EM GCC4.7 ym2=REAL(yPAM) !EM GCC4.7 zm2=REAL(zPAM) !EM GCC4.7 * --------------------------------------------------- * both couples must have a y-cluster * (condition necessary when in RECOVER_SINGLETS mode) * --------------------------------------------------- if(icy1.eq.0.or.icy2.eq.0)goto 111 if(cl_used(icy1).ne.0)goto 111 if(cl_used(icy2).ne.0)goto 111 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - * track parameters on Y VIEW * (2 couples needed) * - - - - - - - - - - - - - - - - - - - - - - - - - - - - if(ndblt.eq.ndblt_max)then if(verbose.eq.1)print*, $ '** warning ** number of identified '// $ 'doublets exceeds vector dimention ' $ ,'( ',ndblt_max,' )' c good2=.false. c goto 880 !fill ntp and go to next event do iv=1,12 c mask_view(iv) = 3 mask_view(iv) = mask_view(iv)+ 2**2 enddo iflag=1 return endif ccc print*,' ',icp1,icp2 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)*(REAL(zini)-zm1)+ym1! EM GCC4.7 zm1, ym1 and alfayz1/2 are REAL **** -----------------------------------------------**** **** reject non phisical couples **** **** -----------------------------------------------**** if(SECOND_SEARCH)goto 111 if( $ abs(alfayz2(ndblt)).gt.alfyz2_max $ .or. $ abs(alfayz1(ndblt)).gt.alfyz1_max $ )ndblt = ndblt-1 111 continue * - - - - - - - - - - - - - - - - - - - - - - - - - - - - * track parameters on Y VIEW - end * - - - - - - - - - - - - - - - - - - - - - - - - - - - - if(icx1.ne.0)then if(cl_used(icx1).ne.0)goto 31 endif if(icx2.ne.0)then if(cl_used(icx2).ne.0)goto 31 endif if(ip2.eq.nplanes)goto 31 !no possible combination with 3 couples do ip3=(ip2+1),nplanes !loop on planes - COPPIA 3 c$$$ print*,'(3) ip ',ip3 c$$$ $ ,mask_view(nviewx(ip3)) c$$$ $ ,mask_view(nviewy(ip3)) if( mask_view(nviewx(ip3)).ne.0 .or. $ mask_view(nviewy(ip3)).ne.0 )goto 30 !skip plane 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$$$ print*,'(3) ip ',ip3,' icp ',icp3 * --------------------------------------------------- * all three couples must have a x-cluster * (condition necessary when in RECOVER_SINGLETS mode) * --------------------------------------------------- if( $ icx1.eq.0.or. $ icx2.eq.0.or. $ icx3.eq.0.or. $ .false.)goto 29 if(cl_used(icx1).ne.0)goto 29 if(cl_used(icx2).ne.0)goto 29 if(cl_used(icx3).ne.0)goto 29 c call xyz_PAM c $ (icx3,icy3,is3,'COG2','COG2',0.,0.)!(1) c call xyz_PAM c $ (icx3,icy3,is3,PFAdef,PFAdef,0.,0.) !(1) call xyz_PAM $ (icx3,icy3,is3,PFAdef,PFAdef $ ,0.,0.,0.,0.) xm3=REAL(xPAM) !EM GCC4.7 ym3=REAL(yPAM) !EM GCC4.7 zm3=REAL(zPAM) !EM GCC4.7 * find the circle passing through the three points iflag_t = DEBUG call tricircle(3,xp,zp,angp,resp,chi $ ,xc,zc,radius,iflag_t) * the circle must intersect the reference plane cc if(iflag.ne.0)goto 29 if(iflag_t.ne.0)then * if tricircle fails, evaluate a straight line if(DEBUG.eq.1) $ print*,'TRICIRCLE failure' $ ,' >>> straight line' radius=0. xc=0. yc=0. SZZ=0. SZX=0. SSX=0. SZ=0. S1=0. X0=0. Ax=0. BX=0. DO I=1,3 XX = XP(I) SZZ=SZZ+ZP(I)*ZP(I) SZX=SZX+ZP(I)*XX SSX=SSX+XX SZ=SZ+ZP(I) S1=S1+1. ENDDO DET=SZZ*S1-SZ*SZ AX=(SZX*S1-SZ*SSX)/DET BX=(SZZ*SSX-SZX*SZ)/DET X0 = AX*REAL(ZINI)+BX ! EM GCC4.7 endif if( .not.SECOND_SEARCH.and. $ radius**2.lt.(ZINI-zc)**2)goto 29 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - * track parameters on X VIEW * (3 couples needed) * - - - - - - - - - - - - - - - - - - - - - - - - - - - - if(ntrpt.eq.ntrpt_max)then if(verbose.eq.1)print*, $ '** warning **'// $ ' number of identified '// $ 'triplets exceeds'// $ ' vector dimention ' $ ,'( ',ntrpt_max,' )' c good2=.false. c goto 880 !fill ntp and go to next event do iv=1,nviews c mask_view(iv) = 4 mask_view(iv) = $ mask_view(iv)+ 2**3 enddo iflag=1 return endif ccc print*,' ',icp1,icp2,icp3 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(radius.ne.0.and.xc.lt.0)then *************POSITIVE DEFLECTION alfaxz1(ntrpt) = xc+sqrt(radius**2-(REAL(ZINI)-zc)**2) !EM GCC4.7 alfaxz2(ntrpt) = (REAL(ZINI)-zc)/ $ sqrt(radius**2-(REAL(ZINI)-zc)**2) !EM GCC4.7 alfaxz3(ntrpt) = 1/radius else if(radius.ne.0.and.xc.ge.0)then *************NEGATIVE DEFLECTION alfaxz1(ntrpt) = xc-sqrt(radius**2-(REAL(ZINI)-zc)**2) alfaxz2(ntrpt) = -(REAL(ZINI)-zc)/ $ sqrt(radius**2-(REAL(ZINI)-zc)**2) !EM GCC4.7 alfaxz3(ntrpt) = -1/radius else if(radius.eq.0)then *************straight fit alfaxz1(ntrpt) = X0 alfaxz2(ntrpt) = AX alfaxz3(ntrpt) = 0. endif c$$$ print*,'alfaxz1 ', alfaxz1(ntrpt) c$$$ print*,'alfaxz2 ', alfaxz2(ntrpt) c$$$ print*,'alfaxz3 ', alfaxz3(ntrpt) **** -----------------------------------------------**** **** reject non phisical triplets **** **** -----------------------------------------------**** if(SECOND_SEARCH)goto 29 if( $ abs(alfaxz2(ntrpt)).gt. $ alfxz2_max $ .or. $ abs(alfaxz1(ntrpt)).gt. $ alfxz1_max $ )ntrpt = ntrpt-1 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - * track parameters on X VIEW - end * - - - - - - - - - - - - - - - - - - - - - - - - - - - - 29 continue enddo !end loop on COPPIA 3 enddo !end loop on sensors - COPPIA 3 30 continue enddo !end loop on planes - COPPIA 3 31 continue c 1 enddo !end loop on COPPIA 2 enddo !end loop on COPPIA 2 enddo !end loop on sensors - COPPIA 2 20 continue enddo !end loop on planes - COPPIA 2 c 11 continue continue enddo !end loop on COPPIA1 enddo !end loop on sensors - COPPIA 1 10 continue enddo !end loop on planes - COPPIA 1 if(DEBUG.EQ.1)then print*,'--- doublets ',ndblt print*,'--- triplets ',ntrpt endif c goto 880 !ntp fill return end *************************************************** * * * * * * * * * * * * ************************************************** subroutine doub_to_YZcloud(iflag) include 'commontracker.f' include 'level1.f' include 'common_momanhough.f' c include 'momanhough_init.f' * 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 if(DEBUG.EQ.1)print*,'doub_to_YZcloud:' *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ * 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 nloop=0 90 continue do idb1=1,ndblt !loop (1) on DOUBLETS if(db_used(idb1).eq.1)goto 2228 !db already included in a cloud 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) if(distance.lt.cutdistyz)then 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) endif 1118 continue enddo !end loop (2) on DOUBLETS c 1188 continue 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).and. $ .true.)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(nplused.lt.nplyz_min)goto 2228 !next doublet * ~~~~~~~~~~~~~~~~~ * >>> NEW CLOUD <<< if(nclouds_yz.ge.ncloyz_max)then if(verbose.eq.1)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 do iv=1,nviews c mask_view(iv) = 5 mask_view(iv) = mask_view(iv) + 2**4 enddo 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) enddo npt_tot=npt_tot+npt if(DEBUG.EQ.1)then 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*,'cpcloud_yz ' $ ,(cpcloud_yz(nclouds_yz,icp),icp=1,ncp_tot) print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes) endif * >>> NEW CLOUD <<< * ~~~~~~~~~~~~~~~~~ 2228 continue enddo !end loop (1) on DOUBLETS if(nloop.lt.nstepy)then cutdistyz = cutdistyz+cutystep nloop = nloop+1 goto 90 endif if(DEBUG.EQ.1)then print*,'Y-Z total clouds ',nclouds_yz endif return end *************************************************** * * * * * * * * * * * * ************************************************** subroutine trip_to_XZcloud(iflag) include 'commontracker.f' include 'level1.f' include 'common_momanhough.f' c include 'momanhough_init.f' * 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 if(DEBUG.EQ.1)print*,'trip_to_XZcloud:' *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ * 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 nloop=0 91 continue do itr1=1,ntrpt !loop (1) on TRIPLETS if(tr_used(itr1).eq.1)goto 22288 !already included in a cloud 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) * ------------------------------------------------------------------------ * FORCE INCLUSION OF TRIPLETS COMPOSED BY SAME COUPLES, IGNORING THE IMAGE * ------------------------------------------------------------------------ * (added in august 2007) istrimage=0 if( $ abs(cpxz1(itrref)).eq.abs(cpxz1(itr2)).and. $ abs(cpxz2(itrref)).eq.abs(cpxz2(itr2)).and. $ abs(cpxz3(itrref)).eq.abs(cpxz3(itr2)).and. $ .true.)istrimage=1 if(distance.lt.cutdistxz.or.istrimage.eq.1)then 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) endif 11188 continue enddo !end loop (2) on TRIPLETS c11888 continue 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 * ------------------------------------------ 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).and. $ .true.)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(nplused.lt.nplxz_min)goto 22288 !next triplet * ~~~~~~~~~~~~~~~~~ * >>> NEW CLOUD <<< if(nclouds_xz.ge.ncloxz_max)then if(verbose.eq.1)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 do iv=1,nviews c mask_view(iv) = 6 mask_view(iv) = mask_view(iv) + 2**5 enddo 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.EQ.1)then 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*,'cpcloud_xz ' $ ,(cpcloud_xz(nclouds_xz,icp),icp=1,ncp_tot) print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes) endif * >>> NEW CLOUD <<< * ~~~~~~~~~~~~~~~~~ 22288 continue enddo !end loop (1) on DOUBLETS if(nloop.lt.nstepx)then cutdistxz=cutdistxz+cutxstep nloop=nloop+1 goto 91 endif if(DEBUG.EQ.1)then print*,'X-Z total clouds ',nclouds_xz endif return end *************************************************** * * * * * * * * * * * * ************************************************** subroutine clouds_to_ctrack(iflag) include 'commontracker.f' include 'level1.f' include 'common_momanhough.f' include 'common_xyzPAM.f' include 'common_mini_2.f' include 'common_mech.f' * 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,2*ncouplemax) integer ncp_match(nplanes) * ----------------------------------------------------------- integer hit_plane(nplanes) * ----------------------------------------------------------- * variables for track fitting double precision AL_INI(5) * ----------------------------------------------------------- if(DEBUG.EQ.1)print*,'clouds_to_ctrack:' ntracks=0 !counter of track candidates do iyz=1,nclouds_yz !loop on YZ clouds do ixz=1,nclouds_xz !loop on XZ clouds * -------------------------------------------------- * 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 !n.matching couples (REAL couples, not SINGLETS) ncp_match(ip)=0 !n.matching couples per plane do icpp=1,ncouplemax cp_match(ip,icpp)=0 !init couple list enddo enddo ncp_ok=0 !count n.matching-couples ncpx_ok=0 !count n.matching-couples with x cluster ncpy_ok=0 !count n.matching-couples with y cluster do icp=1,ncp_tot !loop over couples if(.not.RECOVER_SINGLETS)then * ------------------------------------------------------ * if NOT in RECOVER_SINGLETS mode, take the intersection * between xz yz clouds * ------------------------------------------------------ cpintersec(icp)=min( $ cpcloud_yz(iyz,icp), $ cpcloud_xz(ixz,icp)) * cpintersec is >0 if yz and xz clouds contain the same image of couple icp * ------------------------------------------------------ * discard the couple if the sensor is in conflict * ------------------------------------------------------ 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 else * ------------------------------------------------------ * if RECOVER_SINGLETS take the union * (otherwise the fake couples formed by singlets would be * discarded) * ------------------------------------------------------ cpintersec(icp)=max( $ cpcloud_yz(iyz,icp), $ cpcloud_xz(ixz,icp)) c$$$ if(cpcloud_yz(iyz,icp).gt.0) c$$$ $ cpintersec(icp)=cpcloud_yz(iyz,icp) * cpintersec is >0 if either yz or xz cloud contains the couple icp endif c$$$ print*,icp,ip_cp(icp),' -- ',cpintersec(icp) if(cpintersec(icp).ne.0)then ip=ip_cp(icp) hit_plane(ip)=1 c$$$ if(clx(ip,icp).gt.0.and.cly(ip,icp).gt.0) c$$$ $ ncp_ok=ncp_ok+1 c$$$ if(clx(ip,icp).gt.0.and.cly(ip,icp).eq.0) c$$$ $ ncpx_ok=ncpx_ok+1 c$$$ if(clx(ip,icp).eq.0.and.cly(ip,icp).gt.0) c$$$ $ ncpy_ok=ncpy_ok+1 if( cpcloud_yz(iyz,icp).gt.0.and. $ cpcloud_xz(ixz,icp).gt.0) $ ncp_ok=ncp_ok+1 if( cpcloud_yz(iyz,icp).gt.0.and. $ cpcloud_xz(ixz,icp).eq.0) $ ncpy_ok=ncpy_ok+1 if( cpcloud_yz(iyz,icp).eq.0.and. $ cpcloud_xz(ixz,icp).gt.0) $ ncpx_ok=ncpx_ok+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.3)goto 888 !next combination ccc if(nplused.lt.nplxz_min)goto 888 !next combination ccc if(nplused.lt.nplyz_min)goto 888 !next combination * ----------------------------------------------------------- * if in RECOVER_SINGLET mode, the two clouds must have * at least ONE intersecting real couple * ----------------------------------------------------------- if(ncp_ok.lt.1)goto 888 !next combination if(DEBUG.EQ.1)then print*,'////////////////////////////' print*,'Cloud combination (Y,X): ',iyz,ixz print*,' db ',ptcloud_yz(iyz) print*,' tr ',ptcloud_xz(ixz) print*,' -----> # matching couples ',ncp_ok print*,' -----> # fake couples (X)',ncpx_ok print*,' -----> # fake couples (Y)',ncpy_ok do icp=1,ncp_tot print*,'cp ',icp,' >' $ ,' x ',cpcloud_xz(ixz,icp) $ ,' y ',cpcloud_yz(iyz,icp) $ ,' ==> ',cpintersec(icp) enddo endif if(DEBUG.EQ.1)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 if(DEBUG.eq.1) $ print*,'combination: ' $ ,cp_match(6,icp1) $ ,cp_match(5,icp2) $ ,cp_match(4,icp3) $ ,cp_match(3,icp4) $ ,cp_match(2,icp5) $ ,cp_match(1,icp6) * --------------------------------------- * check if this group of couples has been * already fitted * --------------------------------------- do ica=1,ntracks isthesame=1 do ip=1,NPLANES if(hit_plane(ip).ne.0)then if( CP_STORE(nplanes-ip+1,ica) $ .ne. $ cp_match(ip,hit_plane(ip)) ) $ isthesame=0 else if( CP_STORE(nplanes-ip+1,ica) $ .ne. $ 0 ) $ isthesame=0 endif enddo if(isthesame.eq.1)then if(DEBUG.eq.1) $ print*,'(already fitted)' goto 666 !jump to next combination endif enddo call track_init !init TRACK common do ip=1,nplanes !loop on planes (bottom to top) 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.) c call xyz_PAM(icx,icy,is, !(1) c $ PFAdef,PFAdef,0.,0.) !(1) call xyz_PAM(icx,icy,is, !(1) $ PFAdef,PFAdef,0.,0.,0.,0.) * ************************* * ----------------------------- if(icx.gt.0.and.icy.gt.0)then 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 if(DEBUG.EQ.1)print*,'(X,Y)' $ ,nplanes-ip+1,xPAM,yPAM else xm_A(nplanes-ip+1) = xPAM_A ym_A(nplanes-ip+1) = yPAM_A xm_B(nplanes-ip+1) = xPAM_B ym_B(nplanes-ip+1) = yPAM_B zm(nplanes-ip+1) $ = (zPAM_A+zPAM_B)/2. resx(nplanes-ip+1) = resxPAM resy(nplanes-ip+1) = resyPAM if(icx.eq.0.and.icy.gt.0)then xgood(nplanes-ip+1)=0. ygood(nplanes-ip+1)=1. resx(nplanes-ip+1) = 1000. if(DEBUG.EQ.1)print*,'( Y)' $ ,nplanes-ip+1,xPAM,yPAM elseif(icx.gt.0.and.icy.eq.0)then xgood(nplanes-ip+1)=1. ygood(nplanes-ip+1)=0. if(DEBUG.EQ.1)print*,'(X )' $ ,nplanes-ip+1,xPAM,yPAM resy(nplanes-ip+1) = 1000. else print*,'both icx=0 and icy=0' $ ,' ==> IMPOSSIBLE!!' endif endif * ----------------------------- endif enddo !end loop on planes * ********************************************************** * ************************** FIT *** FIT *** FIT *** FIT *** * ********************************************************** cccc scommentare se si usa al_ini della nuvola c$$$ do i=1,5 c$$$ AL(i)=AL_INI(i) c$$$ enddo call guess() do i=1,5 AL_INI(i)=AL(i) enddo ifail=0 !error flag in chi^2 computation jstep=0 !number of minimization steps iprint=0 c if(DEBUG.EQ.1)iprint=1 if(DEBUG.EQ.1)iprint=2 call mini2(jstep,ifail,iprint) if(ifail.ne.0) then if(DEBUG.EQ.1)then print *, $ '*** MINIMIZATION FAILURE *** ' $ //'(clouds_to_ctrack)' print*,'initial guess: ' print*,'AL_INI(1) = ',AL_INI(1) print*,'AL_INI(2) = ',AL_INI(2) print*,'AL_INI(3) = ',AL_INI(3) print*,'AL_INI(4) = ',AL_INI(4) print*,'AL_INI(5) = ',AL_INI(5) endif c chi2=-chi2 endif * ********************************************************** * ************************** FIT *** FIT *** FIT *** FIT *** * ********************************************************** if(chi2.le.0.)goto 666 if(chi2.ge.1.e08)goto 666 !OPTIMIZATION if(chi2.ne.chi2)goto 666 !OPTIMIZATION * -------------------------- * STORE candidate TRACK INFO * -------------------------- if(ntracks.eq.NTRACKSMAX)then if(verbose.eq.1)print*, $ '** warning ** number of candidate tracks '// $ ' exceeds vector dimension ' $ ,'( ',NTRACKSMAX,' )' c good2=.false. c goto 880 !fill ntp and go to next event do iv=1,nviews c mask_view(iv) = 7 mask_view(iv) = mask_view(iv) + 2**6 enddo iflag=1 return endif ntracks = ntracks + 1 do ip=1,nplanes !top to bottom 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)) * NB! hit_plane is defined from bottom to top if(hit_plane(ip).ne.0)then CP_STORE(nplanes-ip+1,ntracks)= $ cp_match(ip,hit_plane(ip)) SENSOR_STORE(nplanes-ip+1,ntracks) $ = is_cp(cp_match(ip,hit_plane(ip))) icl= $ clx(ip,icp_cp( $ cp_match(ip,hit_plane(ip) $ ))); if(icl.eq.0) $ icl= $ cly(ip,icp_cp( $ cp_match(ip,hit_plane(ip) $ ))); LADDER_STORE(nplanes-ip+1,ntracks) $ = LADDER(icl); else CP_STORE(nplanes-ip+1,ntracks)=0 SENSOR_STORE(nplanes-ip+1,ntracks)=0 LADDER_STORE(nplanes-ip+1,ntracks)=0 endif BX_STORE(ip,ntracks)=0!I dont need it now BY_STORE(ip,ntracks)=0!I dont need it now CLS_STORE(ip,ntracks)=0 do i=1,5 AL_STORE(i,ntracks)=sngl(AL(i)) enddo enddo RCHI2_STORE(ntracks)=REAL(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 cc return endif if(DEBUG.EQ.1)then print*,'****** TRACK CANDIDATES *****************' print*,'# R. chi2 RIG ndof' do i=1,ntracks ndof=0 !(1) do ii=1,nplanes !(1) ndof=ndof !(1) $ +int(xgood_store(ii,i)) !(1) $ +int(ygood_store(ii,i)) !(1) enddo !(1) print*,i,' --- ',rchi2_store(i),' --- ' $ ,1./abs(AL_STORE(5,i)),' --- ',ndof enddo print*,'*****************************************' endif return end *************************************************** * * * * * * * * * * * * ************************************************** subroutine refine_track(ibest) include 'commontracker.f' include 'level1.f' include 'common_momanhough.f' include 'common_xyzPAM.f' include 'common_mini_2.f' include 'common_mech.f' include 'calib.f' * flag to chose PFA character*10 PFA common/FINALPFA/PFA double precision xmm,rxmm,xmm_A,xmm_B !EM GCC4.7 double precision ymm,rymm,ymm_A,ymm_B !EM GCC4.7 double precision zmm,zmm_A,zmm_B !EM GCC4.7 double precision clincnewc !EM GCC4.7 double precision clincnew !EM GCC4.7 real k(6) DATA k/1.099730,0.418900,0.220939,0.220907,0.418771,1.100674/ real xp,yp,zp real xyzp(3),bxyz(3) equivalence (xp,xyzp(1)),(yp,xyzp(2)),(zp,xyzp(3)) if(DEBUG.EQ.1)print*,'refine_track:' * ================================================= * 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(DEBUG.EQ.1)print*,' ........... plane ',ip,' ........... ' xP=XV_STORE(nplanes-ip+1,ibest) yP=YV_STORE(nplanes-ip+1,ibest) zP=ZV_STORE(nplanes-ip+1,ibest) call gufld(xyzp,bxyz) BX_STORE(nplanes-ip+1,ibest)=bxyz(1) BY_STORE(nplanes-ip+1,ibest)=bxyz(2) c$$$ bxyz(1)=0 c$$$ bxyz(2)=0 c$$$ bxyz(3)=0 * ||||||||||||||||||||||||||||||||||||||||||||||||| * ------------------------------------------------- * If the plane has been already included, it just * computes again the coordinates of the x-y couple * using improved PFAs * ------------------------------------------------- * ||||||||||||||||||||||||||||||||||||||||||||||||| c$$$ if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and. c$$$ $ YGOOD_STORE(nplanes-ip+1,ibest).eq.1. )then c$$$ c$$$ id=CP_STORE(nplanes-ip+1,ibest) c$$$ c$$$ is=is_cp(id) c$$$ icp=icp_cp(id) c$$$ if(ip_cp(id).ne.ip) c$$$ $ print*,'OKKIO!!' c$$$ $ ,'id ',id,is,icp c$$$ $ ,ip_cp(id),ip c$$$ icx=clx(ip,icp) c$$$ icy=cly(ip,icp) c$$$c call xyz_PAM(icx,icy,is, c$$$c $ PFA,PFA, c$$$c $ AXV_STORE(nplanes-ip+1,ibest), c$$$c $ AYV_STORE(nplanes-ip+1,ibest)) c$$$ call xyz_PAM(icx,icy,is, c$$$ $ PFA,PFA, c$$$ $ AXV_STORE(nplanes-ip+1,ibest), c$$$ $ AYV_STORE(nplanes-ip+1,ibest), c$$$ $ bxyz(1), c$$$ $ bxyz(2) c$$$ $ ) c$$$ c$$$ xm(nplanes-ip+1) = xPAM c$$$ ym(nplanes-ip+1) = yPAM c$$$ zm(nplanes-ip+1) = zPAM c$$$ xgood(nplanes-ip+1) = 1 c$$$ ygood(nplanes-ip+1) = 1 c$$$ resx(nplanes-ip+1) = resxPAM c$$$ resy(nplanes-ip+1) = resyPAM c$$$ c$$$ dedxtrk_x(nplanes-ip+1)=sgnl(icx)/mip(VIEW(icx),LADDER(icx)) c$$$ dedxtrk_y(nplanes-ip+1)=sgnl(icy)/mip(VIEW(icy),LADDER(icy)) if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..or. $ 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) c call xyz_PAM(icx,icy,is, c $ PFA,PFA, c $ AXV_STORE(nplanes-ip+1,ibest), c $ AYV_STORE(nplanes-ip+1,ibest)) call xyz_PAM(icx,icy,is, $ PFA,PFA, $ AXV_STORE(nplanes-ip+1,ibest), $ AYV_STORE(nplanes-ip+1,ibest), $ bxyz(1), $ bxyz(2) $ ) if(icx.gt.0.and.icy.gt.0)then xm(nplanes-ip+1) = xPAM ym(nplanes-ip+1) = yPAM zm(nplanes-ip+1) = zPAM xm_A(nplanes-ip+1) = 0. ym_A(nplanes-ip+1) = 0. xm_B(nplanes-ip+1) = 0. ym_B(nplanes-ip+1) = 0. xgood(nplanes-ip+1) = 1 ygood(nplanes-ip+1) = 1 resx(nplanes-ip+1) = resxPAM resy(nplanes-ip+1) = resyPAM dedxtrk_x(nplanes-ip+1)= $ sgnl(icx)/mip(VIEW(icx),LADDER(icx)) dedxtrk_y(nplanes-ip+1)= $ sgnl(icy)/mip(VIEW(icy),LADDER(icy)) else xm(nplanes-ip+1) = 0. ym(nplanes-ip+1) = 0. zm(nplanes-ip+1) = (zPAM_A+zPAM_B)/2. xm_A(nplanes-ip+1) = xPAM_A ym_A(nplanes-ip+1) = yPAM_A xm_B(nplanes-ip+1) = xPAM_B ym_B(nplanes-ip+1) = yPAM_B xgood(nplanes-ip+1) = 0 ygood(nplanes-ip+1) = 0 resx(nplanes-ip+1) = 1000.!resxPAM resy(nplanes-ip+1) = 1000.!resyPAM dedxtrk_x(nplanes-ip+1)= 0 dedxtrk_y(nplanes-ip+1)= 0 if(icx.gt.0)then xgood(nplanes-ip+1) = 1 resx(nplanes-ip+1) = resxPAM dedxtrk_x(nplanes-ip+1)= $ sgnl(icx)/mip(VIEW(icx),LADDER(icx)) elseif(icy.gt.0)then ygood(nplanes-ip+1) = 1 resy(nplanes-ip+1) = resyPAM dedxtrk_y(nplanes-ip+1)= $ sgnl(icy)/mip(VIEW(icy),LADDER(icy)) endif endif * ||||||||||||||||||||||||||||||||||||||||||||||||| * ------------------------------------------------- * 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 CP_STORE(nplanes-ip+1,ibest)=0 !re-init CLS_STORE(nplanes-ip+1,ibest)=0 * -------------------------------------------------------------- * determine which ladder and sensor are intersected by the track 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 SENSOR_STORE(nplanes-ip+1,IBEST)=ist LADDER_STORE(nplanes-ip+1,IBEST)=nldt * -------------------------------------------------------------- if(DEBUG.EQ.1)then print*, $ '------ Plane ',ip,' intersected on LADDER ',nldt $ ,' SENSOR ',ist print*, $ '------ coord: ',XP,YP endif * =========================================== * STEP 1 >>>>>>> try to include a new couple * =========================================== distmin=100000000. 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(icx.eq.0.or.icy.eq.0)goto 1188!if fake couple, jump to next if(LADDER(icx).ne.nldt.or. !If the ladder number does not match c $ cl_used(icx).eq.1.or. !or the X cluster is already used c $ cl_used(icy).eq.1.or. !or the Y cluster is already used $ cl_used(icx).ne.0.or. !or the X cluster is already used $ cl_used(icy).ne.0.or. !or the Y cluster is already used $ .false.)goto 1188 !then jump to next couple. * call xyz_PAM(icx,icy,ist, $ PFA,PFA, $ AXV_STORE(nplanes-ip+1,ibest), $ AYV_STORE(nplanes-ip+1,ibest), $ bxyz(1), $ bxyz(2) $ ) distance = distance_to(XP,YP) c distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI id=id_cp(ip,icp,ist) if(DEBUG.EQ.1) $ print*,'( couple ',id $ ,' ) distance ',distance if(distance.lt.distmin)then xmm = xPAM ymm = yPAM zmm = zPAM rxmm = resxPAM rymm = resyPAM distmin = distance idm = id dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2) dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2) clincnewc=10.*dsqrt(rymm**2+rxmm**2 $ +DBLE(RCHI2_STORE(ibest)*k(ip)*(cov(1,1)+cov(2,2))))! EM GCC4.7 endif 1188 continue enddo !end loop on couples on plane icp if(distmin.le.clincnewc)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 !<<< dedxtrk_x(nplanes-ip+1) = dedxmmx !<<< dedxtrk_y(nplanes-ip+1) = dedxmmy !<<< * ----------------------------------- CP_STORE(nplanes-ip+1,ibest)=idm if(DEBUG.EQ.1)print*,'%%%% included couple ',idm $ ,' (dist.= ',distmin,', cut ',clincnewc,' )' goto 133 !next plane endif * ================================================ * STEP 2 >>>>>>> try to include a single cluster * either from a couple or single * ================================================ 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) if(icx.eq.0.or.icy.eq.0)goto 11882!if fake couple, jump to next 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 ----------------------------------------------- c if(cl_used(icx).eq.1)goto 11881 !if the X cluster is already used if(cl_used(icx).ne.0)goto 11881 !if the X cluster is already used !(3) * !jump to the Y cluster c call xyz_PAM(icx,0,ist, c $ PFA,PFA, c $ AXV_STORE(nplanes-ip+1,ibest),0.) call xyz_PAM(icx,0,ist, $ PFA,PFA, $ AXV_STORE(nplanes-ip+1,ibest),0., $ bxyz(1), $ bxyz(2) $ ) distance = distance_to(XP,YP) c distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI if(DEBUG.EQ.1) $ print*,'( cl-X ',icx $ ,' in cp ',id,' ) 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 = sgnl(icx) !(1) dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2) dedxmmy = 0. !(1) endif 11881 continue *----- try cluster y ----------------------------------------------- c if(cl_used(icy).eq.1)goto 11882 !if the Y cluster is already used if(cl_used(icy).ne.0)goto 11882 !if the Y cluster is already used !(3) * !jump to the next couple c call xyz_PAM(0,icy,ist, c $ PFA,PFA, c $ 0.,AYV_STORE(nplanes-ip+1,ibest)) call xyz_PAM(0,icy,ist, $ PFA,PFA, $ 0.,AYV_STORE(nplanes-ip+1,ibest), $ bxyz(1), $ bxyz(2) $ ) distance = distance_to(XP,YP) c distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI if(DEBUG.EQ.1) $ print*,'( cl-Y ',icy $ ,' in cp ',id,' ) 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 = sgnl(icy) !(1) dedxmmx = 0. !(1) dedxmmy = sgnl(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).ne.0.or. !if the cluster is already used !(3) $ 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, $ PFA,PFA, $ AXV_STORE(nplanes-ip+1,ibest),0., $ bxyz(1), $ bxyz(2) $ ) else !<---- Y view call xyz_PAM(0,icl,ist, $ PFA,PFA, $ 0.,AYV_STORE(nplanes-ip+1,ibest), $ bxyz(1), $ bxyz(2) $ ) endif distance = distance_to(XP,YP) c distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI if(DEBUG.EQ.1) $ print*,'( cl-s ',icl $ ,' ) 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 if(mod(VIEW(icl),2).eq.0)then !<---- X view dedxmmx = sgnl(icl)/mip(VIEW(icl),LADDER(icl)) dedxmmy = 0. else !<---- Y view dedxmmx = 0. dedxmmy = sgnl(icl)/mip(VIEW(icl),LADDER(icl)) endif endif 18882 continue enddo !end loop on single clusters if(iclm.ne.0)then if(mod(VIEW(iclm),2).eq.0)then clincnew= $ 20.* !EM GCC4.7 $ dsqrt(rxmm**2 + $ DBLE(RCHI2_STORE(ibest)*k(ip))*cov(1,1)) else if(mod(VIEW(iclm),2).ne.0)then clincnew= $ 10.* !EM GCC4.7 $ dsqrt(rymm**2 + $ DBLE(RCHI2_STORE(ibest)*k(ip))*cov(2,2)) endif if(distmin.le.clincnew)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.EQ.1) $ print*,'%%%% included X-cl ',iclm $ ,'( chi^2, ',RCHI2_STORE(ibest) $ ,', dist.= ',distmin $ ,', cut ',clincnew,' )' else YGOOD(nplanes-ip+1)=1. resy(nplanes-ip+1)=rymm if(DEBUG.EQ.1) $ print*,'%%%% included Y-cl ',iclm $ ,'( chi^2, ',RCHI2_STORE(ibest) $ ,', dist.= ', distmin $ ,', cut ',clincnew,' )' 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. dedxtrk_x(nplanes-ip+1) = dedxmmx !<<< dedxtrk_y(nplanes-ip+1) = dedxmmy !<<< * ---------------------------- endif endif endif 133 continue enddo !end loop on planes return end *************************************************** * * * * * * * * * * * * ************************************************** * * **************************************************** subroutine init_level2 include 'commontracker.f' include 'level1.f' include 'common_momanhough.f' include 'level2.f' * --------------------------------- * variables initialized from level1 * --------------------------------- do i=1,nviews good2(i)=good1(i) do j=1,nva1_view vkflag(i,j)=1 if(cnnev(i,j).le.0)then vkflag(i,j)=cnnev(i,j) endif enddo enddo * ---------------- * level2 variables * ---------------- NTRK = 0 do it=1,NTRKMAX IMAGE(IT)=0 CHI2_nt(IT) = -100000. 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 TAILX_nt(IP,IT) = 0 TAILY_nt(IP,IT) = 0 XBAD(IP,IT) = 0 YBAD(IP,IT) = 0 XGOOD_nt(IP,IT) = 0 YGOOD_nt(IP,IT) = 0 LS(IP,IT) = 0 DEDX_X(IP,IT) = 0 DEDX_Y(IP,IT) = 0 CLTRX(IP,IT) = 0 CLTRY(IP,IT) = 0 multmaxx(ip,it) = 0 seedx(ip,it) = 0 xpu(ip,it) = 0 multmaxy(ip,it) = 0 seedy(ip,it) = 0 ypu(ip,it) = 0 enddo do ipa=1,5 AL_nt(IPA,IT) = 0 do ipaa=1,5 coval(ipa,ipaa,IT)=0 enddo enddo enddo nclsx=0 nclsy=0 do ip=1,NSINGMAX planex(ip)=0 xs(1,ip)=0 xs(2,ip)=0 sgnlxs(ip)=0 planey(ip)=0 ys(1,ip)=0 ys(2,ip)=0 sgnlys(ip)=0 sxbad(ip)=0 sybad(ip)=0 multmaxsx(ip)=0 multmaxsy(ip)=0 enddo end ************************************************************ * * * * * * * ************************************************************ subroutine init_hough include 'commontracker.f' include 'level1.f' include 'common_momanhough.f' include 'common_hough.f' include 'level2.f' 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 alfayz1_nt(idb)=0 alfayz2_nt(idb)=0 enddo do itr=1,ntrpt_max_nt tr_cloud_nt(itr)=0 alfaxz1_nt(itr)=0 alfaxz2_nt(itr)=0 alfaxz3_nt(itr)=0 enddo do idb=1,ncloyz_max ptcloud_yz_nt(idb)=0 alfayz1_av_nt(idb)=0 alfayz2_av_nt(idb)=0 enddo do itr=1,ncloxz_max ptcloud_xz_nt(itr)=0 alfaxz1_av_nt(itr)=0 alfaxz2_av_nt(itr)=0 alfaxz3_av_nt(itr)=0 enddo ntrpt=0 ndblt=0 NCLOUDS_XZ=0 NCLOUDS_YZ=0 do idb=1,ndblt_max db_cloud(idb)=0 cpyz1(idb)=0 cpyz2(idb)=0 alfayz1(idb)=0 alfayz2(idb)=0 enddo do itr=1,ntrpt_max tr_cloud(itr)=0 cpxz1(itr)=0 cpxz2(itr)=0 cpxz3(itr)=0 alfaxz1(itr)=0 alfaxz2(itr)=0 alfaxz3(itr)=0 enddo do idb=1,ncloyz_max ptcloud_yz(idb)=0 alfayz1_av(idb)=0 alfayz2_av(idb)=0 do idbb=1,ncouplemaxtot cpcloud_yz(idb,idbb)=0 enddo enddo do itr=1,ncloxz_max ptcloud_xz(itr)=0 alfaxz1_av(itr)=0 alfaxz2_av(itr)=0 alfaxz3_av(itr)=0 do itrr=1,ncouplemaxtot cpcloud_xz(itr,itrr)=0 enddo 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. * ------------------------------------------------------- include 'commontracker.f' include 'level1.f' include 'common_momanhough.f' include 'level2.f' include 'common_mini_2.f' include 'calib.f' character*10 PFA common/FINALPFA/PFA real sinth,phi,pig, npig ! EM GCC4.7 integer ssensor,sladder pig=acos(-1.) * ------------------------------------- chi2_nt(ntr) = sngl(chi2) nstep_nt(ntr) = nstep * ------------------------------------- phi = REAL(al(4)) sinth = REAL(al(3)) if(sinth.lt.0)then sinth = -sinth phi = phi + pig endif npig = aint(phi/(2*pig)) phi = phi - npig*2*pig if(phi.lt.0) $ phi = phi + 2*pig al(4) = phi al(3) = sinth do i=1,5 al_nt(i,ntr) = sngl(al(i)) do j=1,5 coval(i,j,ntr) = sngl(cov(i,j)) enddo 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)) TAILX_nt(IP,ntr) = 0. TAILY_nt(IP,ntr) = 0. 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)) factor = sqrt( $ tan( acos(-1.) * sngl(axv(ip)) /180. )**2 + $ tan( acos(-1.) * sngl(ayv(ip)) /180. )**2 + $ 1. ) dedx_x(ip,ntr) = sngl(dedxtrk_x(ip)/factor) dedx_y(ip,ntr) = sngl(dedxtrk_y(ip)/factor) ccc print*,ip,'dedx >>> ',dedx_x(ip,ntr),dedx_y(ip,ntr) ax = axv_nt(ip,ntr) ay = ayv_nt(ip,ntr) bfx = BX_STORE(ip,IDCAND) bfy = BY_STORE(ip,IDCAND) c$$$ if(ip.eq.6) ax = -1. * axv_nt(ip,ntr) c$$$ if(ip.eq.6) bfy = -1. * BY_STORE(ip,IDCAND) c$$$ tgtemp = tan(ax*acos(-1.)/180.) + pmuH_h*bfy*0.00001 c$$$ angx = 180.*atan(tgtemp)/acos(-1.) c$$$ tgtemp = tan(ay*acos(-1.)/180.)+pmuH_e*bfx*0.00001 c$$$ angy = 180.*atan(tgtemp)/acos(-1.) angx = effectiveangle(ax,2*ip,bfy) angy = effectiveangle(ay,2*ip-1,bfx) id = CP_STORE(ip,IDCAND) ! couple id icl = CLS_STORE(ip,IDCAND) ssensor = -1 sladder = -1 ssensor = SENSOR_STORE(ip,IDCAND) sladder = LADDER_STORE(ip,IDCAND) if(ip.eq.6.and.ssensor.ne.0)ssensor = 3 - ssensor !notazione paolo x align LS(IP,ntr) = ssensor+10*sladder c if(id.ne.0)then CCCCCC(10 novembre 2009) PATCH X NUCLEI C non ho capito perche', ma durante il ritracciamento dei nuclei C (quando una traccia viene trovata ma non e' buona) c'e' qualche variabile C che non viene reinizializzata correttamente e i cluster esclusi C dal fit risultano ancora inclusi... C cltrx(ip,ntr) = 0 cltry(ip,ntr) = 0 c$$$ if( c$$$ $ xgood_nt(ip,ntr).eq.1.and.ygood_nt(ip,ntr).eq.1 c$$$ $ .and. c$$$ $ id.ne.0)then if(id.ne.0)then !patch 30/12/09 c >>> is a couple cltrx(ip,ntr) = clx(nplanes-ip+1,icp_cp(id)) cltry(ip,ntr) = cly(nplanes-ip+1,icp_cp(id)) if(clx(nplanes-ip+1,icp_cp(id)).gt.0)then cl_used(cltrx(ip,ntr)) = 1 !tag used clusters xbad(ip,ntr)= nbadstrips(4,clx(nplanes-ip+1,icp_cp(id))) if(nsatstrips(clx(nplanes-ip+1,icp_cp(id))).gt.0) $ dedx_x(ip,ntr)=-dedx_x(ip,ntr) multmaxx(ip,ntr) = maxs(cltrx(ip,ntr)) $ +10000*mult(cltrx(ip,ntr)) seedx(ip,ntr) = clsignal(indmax(cltrx(ip,ntr))) $ /clsigma(indmax(cltrx(ip,ntr))) call applypfa(PFA,cltrx(ip,ntr),angx,corr,res) xpu(ip,ntr) = corr endif if(cly(nplanes-ip+1,icp_cp(id)).gt.0)then cl_used(cltry(ip,ntr)) = 1 !tag used clusters ybad(ip,ntr)= nbadstrips(4,cly(nplanes-ip+1,icp_cp(id))) if(nsatstrips(cly(nplanes-ip+1,icp_cp(id))).gt.0) $ dedx_y(ip,ntr)=-dedx_y(ip,ntr) multmaxy(ip,ntr) = maxs(cltry(ip,ntr)) $ +10000*mult(cltry(ip,ntr)) seedy(ip,ntr) = clsignal(indmax(cltry(ip,ntr))) $ /clsigma(indmax(cltry(ip,ntr))) call applypfa(PFA,cltry(ip,ntr),angy,corr,res) ypu(ip,ntr) = corr endif c$$$ elseif(icl.ne.0)then endif !patch 30/12/09 if(icl.ne.0)then !patch 30/12/09 cl_used(icl) = 1 !tag used clusters if(mod(VIEW(icl),2).eq.0)then cltrx(ip,ntr)=icl xbad(ip,ntr) = nbadstrips(4,icl) if(nsatstrips(icl).gt.0)dedx_x(ip,ntr)=-dedx_x(ip,ntr) multmaxx(ip,ntr) = maxs(cltrx(ip,ntr)) $ +10000*mult(cltrx(ip,ntr)) seedx(ip,ntr) = clsignal(indmax(cltrx(ip,ntr))) $ /clsigma(indmax(cltrx(ip,ntr))) call applypfa(PFA,cltrx(ip,ntr),angx,corr,res) xpu(ip,ntr) = corr elseif(mod(VIEW(icl),2).eq.1)then cltry(ip,ntr)=icl ybad(ip,ntr) = nbadstrips(4,icl) if(nsatstrips(icl).gt.0)dedx_y(ip,ntr)=-dedx_y(ip,ntr) multmaxy(ip,ntr) = maxs(cltry(ip,ntr)) $ +10000*mult(cltry(ip,ntr)) seedy(ip,ntr) = clsignal(indmax(cltry(ip,ntr))) $ /clsigma(indmax(cltry(ip,ntr))) call applypfa(PFA,cltry(ip,ntr),angy,corr,res) ypu(ip,ntr) = corr endif endif enddo if(DEBUG.eq.1)then print*,'> STORING TRACK ',ntr print*,'clusters: ' do ip=1,6 print*,'> ',ip,' -- ',cltrx(ip,ntr),cltry(ip,ntr) enddo print*,'dedx: ' do ip=1,6 print*,'> ',ip,' -- ',dedx_x(ip,ntr),dedx_y(ip,ntr) enddo endif end subroutine fill_level2_siglets * ------------------------------------------------------- * 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 'commontracker.f' include 'calib.f' include 'level1.f' include 'common_momanhough.f' include 'level2.f' include 'common_xyzPAM.f' * count #cluster per plane not associated to any track nclsx = 0 nclsy = 0 do iv = 1,nviews c if( mask_view(iv).ne.0 )good2(iv) = 20+mask_view(iv) good2(iv) = good2(iv) + mask_view(iv)*2**8 enddo if(DEBUG.eq.1)then print*,'> STORING SINGLETS ' endif do icl=1,nclstr1 ip=nplanes-npl(VIEW(icl))+1 if(cl_used(icl).eq.0)then !cluster not included in any track if(mod(VIEW(icl),2).eq.0)then !=== X views nclsx = nclsx + 1 planex(nclsx) = ip sgnlxs(nclsx) = sgnl(icl)/mip(VIEW(icl),LADDER(icl)) if(nsatstrips(icl).gt.0)sgnlxs(nclsx)=-sgnlxs(nclsx) clsx(nclsx) = icl sxbad(nclsx) = nbadstrips(1,icl) multmaxsx(nclsx) = maxs(icl)+10000*mult(icl) do is=1,2 c call xyz_PAM(icl,0,is,'COG1',' ',0.,0.) c call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.) call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.,0.,0.) xs(is,nclsx) = REAL((xPAM_A+xPAM_B)/2.) ! EM GCC4.7 enddo else !=== Y views nclsy = nclsy + 1 planey(nclsy) = ip sgnlys(nclsy) = sgnl(icl)/mip(VIEW(icl),LADDER(icl)) if(nsatstrips(icl).gt.0)sgnlys(nclsy)=-sgnlys(nclsy) clsy(nclsy) = icl sybad(nclsy) = nbadstrips(1,icl) multmaxsy(nclsy) = maxs(icl)+10000*mult(icl) do is=1,2 c call xyz_PAM(0,icl,is,' ','COG1',0.,0.) c call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.) call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.,0.,0.) ys(is,nclsy) = REAL((yPAM_A+yPAM_B)/2.) ! EM GCC4.7 enddo endif endif ***** LO METTO QUI PERCHE` NON SO DOVE METTERLO whichtrack(icl) = cl_used(icl) * -------------------------------------------------- * per non perdere la testa... * whichtrack(icl) e` una variabile del common level1 * che serve solo per sapere quali cluster sono stati * associati ad una traccia, e permettere di salvare * solo questi nell'albero di uscita * -------------------------------------------------- enddo end *************************************************** * * * * * * * * * * * * ************************************************** subroutine fill_hough * ------------------------------------------------------- * This routine fills the variables related to the hough * transform, for the debig n-tuple * ------------------------------------------------------- include 'commontracker.f' include 'level1.f' include 'common_momanhough.f' include 'common_hough.f' include 'level2.f' 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 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 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 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,min(ndblt_max_nt,nnn) db_cloud_nt(ipt)=db_cloud(ipt) enddo 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,min(ntrpt_max_nt,nnn) tr_cloud_nt(ipt)=tr_cloud(ipt) enddo endif endif end