--- DarthVader/TrackerLevel2/src/F77/analysissubroutines.f 2006/11/02 11:19:47 1.10
+++ DarthVader/TrackerLevel2/src/F77/analysissubroutines.f 2014/01/16 15:29:50 1.45
@@ -18,10 +18,40 @@
include 'common_xyzPAM.f'
include 'common_mini_2.f'
include 'calib.f'
-c include 'level1.f'
include 'level2.f'
-c include 'momanhough_init.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
@@ -42,12 +72,11 @@
*-------------------------------------------------------------------------------
*-------------------------------------------------------------------------------
-c iflag=0
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
@@ -76,11 +105,12 @@
*-------------------------------------------------------------------------------
*-------------------------------------------------------------------------------
-c iflag=0
+
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
*-------------------------------------------------------------------------------
@@ -131,37 +161,65 @@
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
- goto 878
- endif
+ 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(planehit.eq.3) goto 881
+ 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
- 881 continue
-* if there is at least three planes on the Y view decreases cuts on X view
- if(nclouds_xz.eq.0.and.nclouds_yz.gt.0.and.
- $ nplxz_min.ne.y_min_start)then
- nptxz_min=x_min_step
- nplxz_min=x_min_start-x_min_step
- goto 879
- endif
+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
@@ -183,6 +241,8 @@
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)
@@ -209,15 +269,20 @@
*
*-------------------------------------------------------------------------------
*-------------------------------------------------------------------------------
- ntrk=0 !counter of identified physical tracks
+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
-11111 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
@@ -233,24 +298,35 @@
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 !(1)
+ ndofbest=0
do i=1,ntracks
- if(RCHI2_STORE(i).lt.rchi2best.and.
- $ RCHI2_STORE(i).gt.0)then
- 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)
- if(ndof.ge.ndofbest)then !(1)
+ 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 !(1)
- endif !(1)
+ 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.
@@ -269,8 +345,10 @@
iimage=0
endif
if(icand.eq.0)then
- print*,'HAI FATTO UN CASINO!!!!!! icand = ',icand
- $ ,ibest,iimage
+ if(VERBOSE.EQ.1)then
+ print*,'HAI FATTO UN CASINO!!!!!! icand = ',icand
+ $ ,ibest,iimage
+ endif
return
endif
@@ -281,26 +359,35 @@
* **********************************************************
* ************************** 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
- if(DEBUG)iprint=1
+c if(DEBUG.EQ.1)iprint=1
+ if(VERBOSE.EQ.1)iprint=1
+ if(DEBUG.EQ.1)iprint=2
call mini2(jstep,ifail,iprint)
- if(ifail.ne.0) then
- if(DEBUG)then
+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 *** (mini2) '
+ $ '*** MINIMIZATION FAILURE *** (after refinement) '
$ ,iev
endif
- chi2=-chi2
endif
- if(DEBUG)then
+ if(DEBUG.EQ.1)then
print*,'----------------------------- improved track coord'
22222 format(i2,' * ',3f10.4,' --- ',4f10.4,' --- ',2f4.0,2f10.5)
do ip=1,6
@@ -311,7 +398,7 @@
endif
c rchi2=chi2/dble(ndof)
- if(DEBUG)then
+ if(DEBUG.EQ.1)then
print*,' '
print*,'****** SELECTED TRACK *************'
print*,'# R. chi2 RIG'
@@ -327,38 +414,122 @@
* ------------- search if the track has an IMAGE -------------
* ------------- (also this is stored ) -------------
if(FIMAGE)goto 122 !>>> jump! (this is already an image)
-* now search for track-image, by comparing couples IDs
+
+* -----------------------------------------------------
+* 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
- iimage=i
+ ncp=0
do ip=1,nplanes
- if( CP_STORE(nplanes-ip+1,icand).ne.
- $ -1*CP_STORE(nplanes-ip+1,i) )iimage=0
+ 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( iimage.ne.0.and.
-c $ RCHI2_STORE(i).le.CHI2MAX.and.
-c $ RCHI2_STORE(i).gt.0.and.
- $ .true.)then
- if(DEBUG)print*,'Track candidate ',iimage
+
+ if(DEBUG.EQ.1)then
+ print*,'Track candidate ',iimage
$ ,' >>> TRACK IMAGE >>> of'
$ ,ibest
- goto 122 !image track found
endif
- enddo
+
+ endif
+
+
122 continue
-* --- and store the results --------------------------------
- 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 print*,'++++++++++ iimage,fimage,ntrk,image '
-c $ ,iimage,fimage,ntrk,image(ntrk)
+* --- and store the results --------------------------------
if(ntrk.eq.NTRKMAX)then
- if(verbose)
+ if(verbose.eq.1)
$ print*,
$ '** warning ** number of identified '//
$ 'tracks exceeds vector dimension '
@@ -366,220 +537,34 @@
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
-* --- then remove selected clusters (ibest+iimage) from clouds ----
- call clean_XYclouds(ibest,iflag)
- if(iflag.eq.1)then !bad event
- goto 880 !fill ntp and go to next event
- endif
-
-* **********************************************************
-* condition to start a new search
-* **********************************************************
- ixznew=0
- do ixz=1,nclouds_xz
- if(ptcloud_xz(ixz).ge.nptxz_min)ixznew=1
- enddo
- iyznew=0
- do iyz=1,nclouds_yz
- if(ptcloud_yz(iyz).ge.nptyz_min)iyznew=1
- enddo
-
- if(ixznew.ne.0.and.
- $ iyznew.ne.0.and.
- $ rchi2best.le.CHI2MAX.and.
-c $ rchi2best.lt.15..and.
- $ .true.)then
- if(DEBUG)then
- print*,'***** NEW SEARCH ****'
- endif
- goto 11111 !try new search
-
- endif
-* **********************************************
-
-
880 return
end
-
-
-c$$$************************************************************
-c$$$
-c$$$ subroutine readmipparam
-c$$$
-c$$$ include 'commontracker.f'
-c$$$ include 'calib.f'
-c$$$
-c$$$ character*60 fname_param
-c$$$ 201 format('trk-LADDER',i1,'-mip.dat')
-c$$$ do ilad=1,nladders_view
-c$$$ write(fname_param,201)ilad
-c$$$ print *,'Opening file: ',fname_param
-c$$$ open(10,
-c$$$ $ FILE='./bin-aux/'//fname_param(1:LNBLNK(fname_param))
-c$$$ $ ,STATUS='UNKNOWN'
-c$$$ $ ,IOSTAT=iostat
-c$$$ $ )
-c$$$ if(iostat.ne.0)then
-c$$$ print*,'READMIPPARAM: *** Error in opening file ***'
-c$$$ return
-c$$$ endif
-c$$$ do iv=1,nviews
-c$$$ read(10,*
-c$$$ $ ,IOSTAT=iostat
-c$$$ $ )pip,
-c$$$ $ mip(int(pip),ilad)
-c$$$c print*,ilad,iv,pip,mip(int(pip),ilad)
-c$$$ enddo
-c$$$ close(10)
-c$$$ enddo
-c$$$
-c$$$ return
-c$$$ end
-c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
-c$$$ subroutine readchargeparam
-c$$$
-c$$$
-c$$$ include 'commontracker.f'
-c$$$ include 'calib.f'
-c$$$
-c$$$ character*60 fname_param
-c$$$ 201 format('charge-l',i1,'.dat')
-c$$$ do ilad=1,nladders_view
-c$$$ write(fname_param,201)ilad
-c$$$ print *,'Opening file: ',fname_param
-c$$$ open(10,
-c$$$ $ FILE='./bin-aux/'//fname_param(1:LNBLNK(fname_param))
-c$$$ $ ,STATUS='UNKNOWN'
-c$$$ $ ,IOSTAT=iostat
-c$$$ $ )
-c$$$ if(iostat.ne.0)then
-c$$$ print*,'READCHARGEPARAM: *** Error in opening file ***'
-c$$$ return
-c$$$ endif
-c$$$ do ip=1,nplanes
-c$$$ read(10,*
-c$$$ $ ,IOSTAT=iostat
-c$$$ $ )pip,
-c$$$ $ kch(ip,ilad),cch(ip,ilad),sch(ip,ilad)
-c$$$c print*,ilad,ip,pip,kch(ip,ilad),
-c$$$c $ cch(ip,ilad),sch(ip,ilad)
-c$$$ enddo
-c$$$ close(10)
-c$$$ enddo
-c$$$
-c$$$ return
-c$$$ end
-c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
-c$$$ subroutine readetaparam
-c$$$* -----------------------------------------
-c$$$* read eta2,3,4 calibration parameters
-c$$$* and fill variables:
-c$$$*
-c$$$* eta2(netabin,nladders_view,nviews)
-c$$$* eta3(2*netabin,nladders_view,nviews)
-c$$$* eta4(2*netabin,nladders_view,nviews)
-c$$$*
-c$$$ include 'commontracker.f'
-c$$$ include 'calib.f'
-c$$$
-c$$$ character*40 fname_binning
-c$$$ character*40 fname_param
-c$$$c character*120 cmd1
-c$$$c character*120 cmd2
-c$$$
-c$$$
-c$$$******retrieve ANGULAR BINNING info
-c$$$ fname_binning='binning.dat'
-c$$$ print *,'Opening file: ',fname_binning
-c$$$ open(10,
-c$$$ $ FILE='./bin-aux/'//fname_binning(1:LNBLNK(fname_binning))
-c$$$ $ ,STATUS='UNKNOWN'
-c$$$ $ ,IOSTAT=iostat
-c$$$ $ )
-c$$$ if(iostat.ne.0)then
-c$$$ print*,'READETAPARAM: *** Error in opening file ***'
-c$$$ return
-c$$$ endif
-c$$$ print*,'---- ANGULAR BINNING ----'
-c$$$ print*,'Bin - angL - angR'
-c$$$ 101 format(i2,' ',f6.2,' ',f6.2)
-c$$$ do ibin=1,nangmax
-c$$$ read(10,*
-c$$$ $ ,IOSTAT=iostat
-c$$$ $ )xnn,angL(ibin),angR(ibin)
-c$$$ if(iostat.ne.0)goto 1000
-c$$$ write(*,101)int(xnn),angL(ibin),angR(ibin)
-c$$$ enddo
-c$$$ 1000 nangbin=int(xnn)
-c$$$ close(10)
-c$$$ print*,'-------------------------'
-c$$$
-c$$$
-c$$$
-c$$$ do ieta=2,4 !loop on eta 2,3,4
-c$$$******retrieve correction parameters
-c$$$ 200 format(' Opening eta',i1,' files...')
-c$$$ write(*,200)ieta
-c$$$
-c$$$ 201 format('eta',i1,'-bin',i1,'-l',i1,'.dat')
-c$$$ 202 format('eta',i1,'-bin',i2,'-l',i1,'.dat')
-c$$$ do iang=1,nangbin
-c$$$ do ilad=1,nladders_view
-c$$$ if(iang.lt.10)write(fname_param,201)ieta,iang,ilad
-c$$$ if(iang.ge.10)write(fname_param,202)ieta,iang,ilad
-c$$$c print *,'Opening file: ',fname_param
-c$$$ open(10,
-c$$$ $ FILE='./bin-aux/'//fname_param(1:LNBLNK(fname_param))
-c$$$ $ ,STATUS='UNKNOWN'
-c$$$ $ ,IOSTAT=iostat
-c$$$ $ )
-c$$$ if(iostat.ne.0)then
-c$$$ print*,'READETAPARAM: *** Error in opening file ***'
-c$$$ return
-c$$$ endif
-c$$$ do ival=1,netavalmax
-c$$$ if(ieta.eq.2)read(10,*
-c$$$ $ ,IOSTAT=iostat
-c$$$ $ )
-c$$$ $ eta2(ival,iang),
-c$$$ $ (feta2(ival,iv,ilad,iang),iv=1,nviews)
-c$$$ if(ieta.eq.3)read(10,*
-c$$$ $ ,IOSTAT=iostat
-c$$$ $ )
-c$$$ $ eta3(ival,iang),
-c$$$ $ (feta3(ival,iv,ilad,iang),iv=1,nviews)
-c$$$ if(ieta.eq.4)read(10,*
-c$$$ $ ,IOSTAT=iostat
-c$$$ $ )
-c$$$ $ eta4(ival,iang),
-c$$$ $ (feta4(ival,iv,ilad,iang),iv=1,nviews)
-c$$$ if(iostat.ne.0)then
-c$$$ netaval=ival-1
-c$$$c$$$ if(eta2(1,iang).ne.-eta2(netaval,iang))
-c$$$c$$$ $ print*,'**** ERROR on parameters !!! ****'
-c$$$ goto 2000
-c$$$ endif
-c$$$ enddo
-c$$$ 2000 close(10)
-c$$$* print*,'... done'
-c$$$ enddo
-c$$$ enddo
-c$$$
-c$$$ enddo !end loop on eta 2,3,4
-c$$$
-c$$$
-c$$$ return
-c$$$ end
-c$$$
-
************************************************************
************************************************************
@@ -604,6 +589,8 @@
* 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.
@@ -642,179 +629,144 @@
*
*
- subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,angx,angy)
+ subroutine xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
-c*****************************************************
-c 07/10/2005 modified by elena vannuccini --> (1)
-c 01/02/2006 modified by elena vannuccini --> (2)
-c 02/02/2006 modified by Elena Vannuccini --> (3)
-c (implemented new p.f.a.)
-c 03/02/2006 modified by Elena Vannuccini --> (4)
-c (implemented variable resolution)
-c*****************************************************
include 'commontracker.f'
include 'level1.f'
include 'calib.f'
-c include 'level1.f'
include 'common_align.f'
include 'common_mech.f'
include 'common_xyzPAM.f'
- include 'common_resxy.f'
-
-c logical DEBUG
-c common/dbg/DEBUG
integer icx,icy !X-Y cluster ID
integer sensor
integer viewx,viewy
character*4 PFAx,PFAy !PFA to be used
- real angx,angy !X-Y angle
+ real 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
-c double precision xi,yi,zi
-c double precision xi_A,yi_A,zi_A
-c double precision xi_B,yi_B,zi_B
parameter (ndivx=30)
+
+
resxPAM = 0
resyPAM = 0
- xPAM = 0.
- yPAM = 0.
- zPAM = 0.
- xPAM_A = 0.
- yPAM_A = 0.
- zPAM_A = 0.
- xPAM_B = 0.
- yPAM_B = 0.
- zPAM_B = 0.
+ 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))
- resxPAM = RESXAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!!
-
- stripx = float(MAXS(icx))
- if(PFAx.eq.'COG1')then !(1)
- stripx = stripx !(1)
- resxPAM = resxPAM !(1)
- elseif(PFAx.eq.'COG2')then
- stripx = stripx + cog(2,icx)
- resxPAM = resxPAM*fbad_cog(2,icx)
- elseif(PFAx.eq.'ETA2')then
-c cog2 = cog(2,icx)
-c etacorr = pfaeta2(cog2,viewx,nldx,angx)
-c stripx = stripx + etacorr
- stripx = stripx + pfaeta2(icx,angx) !(3)
- resxPAM = risx_eta2(angx) ! (4)
- if(DEBUG.and.fbad_cog(2,icx).ne.1)
- $ print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
- resxPAM = resxPAM*fbad_cog(2,icx)
- elseif(PFAx.eq.'ETA3')then !(3)
- stripx = stripx + pfaeta3(icx,angx) !(3)
- resxPAM = risx_eta3(angx) ! (4)
- if(DEBUG.and.fbad_cog(3,icx).ne.1) !(3)
- $ print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)!(3)
- resxPAM = resxPAM*fbad_cog(3,icx) !(3)
- elseif(PFAx.eq.'ETA4')then !(3)
- stripx = stripx + pfaeta4(icx,angx) !(3)
- resxPAM = risx_eta4(angx) ! (4)
- if(DEBUG.and.fbad_cog(4,icx).ne.1) !(3)
- $ print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)!(3)
- resxPAM = resxPAM*fbad_cog(4,icx) !(3)
- elseif(PFAx.eq.'ETA')then !(3)
- stripx = stripx + pfaeta(icx,angx) !(3)
- resxPAM = ris_eta(icx,angx) ! (4)
- if(DEBUG.and.fbad_cog(2,icx).ne.1) !(3)
- $ print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)!(3)
-c resxPAM = resxPAM*fbad_cog(2,icx) !(3)TEMPORANEO
- resxPAM = resxPAM*fbad_eta(icx,angx) !(3)(4)
- elseif(PFAx.eq.'COG')then !(2)
- stripx = stripx + cog(0,icx) !(2)
- resxPAM = risx_cog(angx) ! (4)
- resxPAM = resxPAM*fbad_cog(0,icx)!(2)
- else
- print*,'*** Non valid p.f.a. (x) --> ',PFAx
+
+ 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
-c if(icy.eq.0.and.icx.ne.0)
-c $ print*,PFAx,icx,angx,stripx,resxPAM,'***'
-
+
* -----------------
* CLUSTER Y
* -----------------
if(icy.ne.0)then
+
viewy = VIEW(icy)
nldy = nld(MAXS(icy),VIEW(icy))
nply = npl(VIEW(icy))
- resyPAM = RESYAV !!!!!!!TEMPORANEO!!!!!!!!!!!!!!!!
+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
- print*,'xyz_PAM ***ERROR*** invalid cluster couple!!! '
- $ ,icx,icy
+ 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)
- stripy = float(MAXS(icy))
- if(PFAy.eq.'COG1')then !(1)
- stripy = stripy !(1)
- resyPAM = resyPAM !(1)
- elseif(PFAy.eq.'COG2')then
- stripy = stripy + cog(2,icy)
- resyPAM = resyPAM*fbad_cog(2,icy)
- elseif(PFAy.eq.'ETA2')then
-c cog2 = cog(2,icy)
-c etacorr = pfaeta2(cog2,viewy,nldy,angy)
-c stripy = stripy + etacorr
- stripy = stripy + pfaeta2(icy,angy) !(3)
- resyPAM = risy_eta2(angy) ! (4)
- resyPAM = resyPAM*fbad_cog(2,icy)
- if(DEBUG.and.fbad_cog(2,icy).ne.1)
- $ print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
- elseif(PFAy.eq.'ETA3')then !(3)
- stripy = stripy + pfaeta3(icy,angy) !(3)
- resyPAM = resyPAM*fbad_cog(3,icy) !(3)
- if(DEBUG.and.fbad_cog(3,icy).ne.1) !(3)
- $ print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)!(3)
- elseif(PFAy.eq.'ETA4')then !(3)
- stripy = stripy + pfaeta4(icy,angy) !(3)
- resyPAM = resyPAM*fbad_cog(4,icy) !(3)
- if(DEBUG.and.fbad_cog(4,icy).ne.1) !(3)
- $ print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)!(3)
- elseif(PFAy.eq.'ETA')then !(3)
- stripy = stripy + pfaeta(icy,angy) !(3)
- resyPAM = ris_eta(icy,angy) ! (4)
-c resyPAM = resyPAM*fbad_cog(2,icy) !(3)TEMPORANEO
- resyPAM = resyPAM*fbad_eta(icy,angy) ! (4)
- if(DEBUG.and.fbad_cog(2,icy).ne.1) !(3)
- $ print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)!(3)
- elseif(PFAy.eq.'COG')then
- stripy = stripy + cog(0,icy)
- resyPAM = risy_cog(angy) ! (4)
-c resyPAM = ris_eta(icy,angy) ! (4)
- resyPAM = resyPAM*fbad_cog(0,icy)
- else
- print*,'*** Non valid p.f.a. (x) --> ',PFAx
- endif
+ call applypfa(PFAy,icy,angy,corr,res)
+ stripy = stripy + corr
+ resyPAM = res
+ 20 continue
endif
-
+
c===========================================================
C COUPLE
C===========================================================
@@ -825,14 +777,15 @@
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...
- print*,'xyz_PAM (couple):',
- $ ' WARNING: false X strip: strip ',stripx
+ if(DEBUG.EQ.1) then
+ print*,'xyz_PAM (couple):',
+ $ ' WARNING: false X strip: strip ',stripx
+ endif
endif
- xi = acoordsi(stripx,viewx)
- yi = acoordsi(stripy,viewy)
- zi = 0.
+ xi = dcoordsi(stripx,viewx)
+ yi = dcoordsi(stripy,viewy)
+ zi = 0.D0
-
c------------------------------------------------------------------------
c (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame
c------------------------------------------------------------------------
@@ -866,13 +819,13 @@
yPAM = dcoord(yrt,viewy,nldy,sensor) / 1.d4
zPAM = ( zrt + z_mech_sensor(nplx,nldx,sensor)*1000. ) / 1.d4
- xPAM_A = 0.
- yPAM_A = 0.
- zPAM_A = 0.
-
- xPAM_B = 0.
- yPAM_B = 0.
- zPAM_B = 0.
+ 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.
@@ -892,7 +845,9 @@
nldx = nldy
viewx = viewy + 1
- yi = acoordsi(stripy,viewy)
+ 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
@@ -902,8 +857,6 @@
yi_B = yi
zi_B = 0.
-c print*,'Y-cl ',icy,stripy,' --> ',yi
-c print*,xi_A,' <--> ',xi_B
elseif(icx.ne.0)then
c===========================================================
@@ -914,14 +867,17 @@
nldy = nldx
viewy = viewx - 1
-c print*,'X-singlet ',icx,nplx,nldx,viewx,stripx
-c if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...
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...
- print*,'xyz_PAM (X-singlet):',
- $ ' WARNING: false X strip: strip ',stripx
+ if(DEBUG.EQ.1) then
+ print*,'xyz_PAM (X-singlet):',
+ $ ' WARNING: false X strip: strip ',stripx
+ endif
endif
- xi = acoordsi(stripx,viewx)
+
+ 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
@@ -937,14 +893,13 @@
yi_B = yi
endif
-c print*,'X-cl ',icx,stripx,' --> ',xi
-c print*,yi_A,' <--> ',yi_B
else
-
- print *,'routine xyz_PAM ---> not properly used !!!'
- print *,'icx = ',icx
- print *,'icy = ',icy
+ if(DEBUG.EQ.1) then
+ print *,'routine xyz_PAM ---> not properly used !!!'
+ print *,'icx = ',icx
+ print *,'icy = ',icy
+ endif
goto 100
endif
@@ -983,6 +938,24 @@
$ + 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
@@ -993,9 +966,12 @@
c in PAMELA reference system
c------------------------------------------------------------------------
- xPAM = 0.
- yPAM = 0.
- zPAM = 0.
+ 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
@@ -1006,19 +982,191 @@
zPAM_B = ( zrt_B + z_mech_sensor(nplx,nldx,sensor)*1000.)/ 1.d4
-c print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')'
else
-
- print *,'routine xyz_PAM ---> not properly used !!!'
- print *,'icx = ',icx
- print *,'icy = ',icy
-
+ 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
********************************************************************************
********************************************************************************
@@ -1040,7 +1188,7 @@
*
********************************************************************************
- real function distance_to(XPP,YPP)
+ real function distance_to(rXPP,rYPP)
include 'common_xyzPAM.f'
@@ -1049,14 +1197,19 @@
* ( 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 (
- + xPAM.eq.0.and.
- + yPAM.eq.0.and.
- + zPAM.eq.0.and.
+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.
@@ -1094,19 +1247,17 @@
endif
distance=
- $ ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2
+ $ ((xmi-XPP)**2+(ymi-YPP)**2)!QUIQUI
+cc $ ((xmi-XPP)**2+(ymi-YPP)**2)/RE**2
distance=dsqrt(distance)
-c$$$ print*,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b
-c$$$ $ ,' --- distance_to --- ',xpp,ypp
-c$$$ print*,' resolution ',re
* ----------------------
elseif(
- + xPAM.ne.0.and.
- + yPAM.ne.0.and.
- + zPAM.ne.0.and.
+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.
@@ -1119,22 +1270,17 @@
* ----------------------
distance=
- $ ((xPAM-XPP)/resxPAM)**2
- $ +
- $ ((yPAM-YPP)/resyPAM)**2
+ $ ((xPAM-XPP))**2 !QUIQUI
+ $ +
+ $ ((yPAM-YPP))**2
+c$$$ $ ((xPAM-XPP)/resxPAM)**2
+c$$$ $ +
+c$$$ $ ((yPAM-YPP)/resyPAM)**2
distance=dsqrt(distance)
-c$$$ print*,xPAM,yPAM,zPAM
-c$$$ $ ,' --- distance_to --- ',xpp,ypp
-c$$$ print*,' resolution ',resxPAM,resyPAM
else
- print*
- $ ,' function distance_to ---> wrong usage!!!'
- print*,' xPAM,yPAM,zPAM ',xPAM,yPAM,zPAM
- print*,' xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b '
- $ ,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b
endif
distance_to = sngl(distance)
@@ -1174,17 +1320,17 @@
data c1/1.,0.,0.,1./
data c2/1.,-1.,-1.,1./
data c3/1.,1.,0.,0./
- real*8 yvvv,xvvv
+ double precision yvvv,xvvv
double precision xi,yi,zi
double precision xrt,yrt,zrt
real AA,BB
- real yvv(4),xvv(4)
+ double precision yvv(4),xvv(4)
* tollerance to consider the track inside the sensitive area
real ptoll
data ptoll/150./ !um
- external nviewx,nviewy,acoordsi,dcoord
+ external nviewx,nviewy,dcoordsi,dcoord
nplpt = nplPAM !plane
viewx = nviewx(nplpt)
@@ -1199,15 +1345,9 @@
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...
-c if((stripx.le.3).or.(stripx.ge.1022)) then !X has 1018 strips...
- print*,'whichsensor: ',
- $ ' WARNING: false X strip: strip ',stripx
- endif
- xi = acoordsi(stripx,viewx)
- yi = acoordsi(stripy,viewy)
- zi = 0.
+ xi = dcoordsi(stripx,viewx)
+ yi = dcoordsi(stripy,viewy)
+ zi = 0.D0
c------------------------------------------------------------------------
c (xrt,yrt,zrt) = rototranslated coordinates in the silicon sensor frame
c------------------------------------------------------------------------
@@ -1232,8 +1372,6 @@
yvv(iv)=sngl(yvvv)
xvv(iv)=sngl(xvvv)
-c print*,'LADDER ',il,' SENSOR ',is,' vertexes >> '
-c $ ,iv,xvv(iv),yvv(iv)
enddo !end loop on sensor vertexes
dtot=0.
@@ -1241,8 +1379,8 @@
iv1=iside
iv2=mod(iside,4)+1
* straight line passing trhough two consecutive vertexes
- AA = (yvv(iv1)-yvv(iv2))/(xvv(iv1)-xvv(iv2))
- BB = yvv(iv1) - AA*xvv(iv1)
+ 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
@@ -1254,8 +1392,8 @@
iv1=iside
iv2=mod(iside,4)+1
* straight line passing trhough two consecutive vertexes
- AA = (xvv(iv1)-xvv(iv2))/(yvv(iv1)-yvv(iv2))
- BB = xvv(iv1) - AA*yvv(iv1)
+ 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
@@ -1358,7 +1496,6 @@
is_cp=0
if(id.lt.0)is_cp=1
if(id.gt.0)is_cp=2
- if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!'
return
end
@@ -1429,274 +1566,6 @@
*************************************************************************
*************************************************************************
*************************************************************************
-c$$$ subroutine book_debug
-c$$$
-c$$$ include 'commontracker.f'
-c$$$ include 'common_momanhough.f'
-c$$$ include 'common_level2debug.f'
-c$$$
-c$$$ character*35 block1,block2,block3!,block4
-c$$$ $ ,block5!,block6
-c$$$
-c$$$* * * * * * * * * * * * * * * * * * * * * * * *
-c$$$* HOUGH TRANSFORM PARAMETERS
-c$$$
-c$$$ call HBOOK2(1003
-c$$$ $ ,'y vs tg thyz'
-c$$$ $ ,300,-1.,1. !x
-c$$$ $ ,3000,-70.,70.,0.) !y
-c$$$
-c$$$ call HBOOK1(1004
-c$$$ $ ,'Dy'
-c$$$ $ ,100,0.,2.,0.)
-c$$$
-c$$$ call HBOOK1(1005
-c$$$ $ ,'D thyz'
-c$$$ $ ,100,0.,.05,0.)
-c$$$
-c$$$
-c$$$
-c$$$* DEBUG ntuple:
-c$$$ call HBNT(ntp_level2+1,'LEVEL2',' ')
-c$$$
-c$$$ call HBNAME(ntp_level2+1,'EVENT',good2_nt,
-c$$$ $ 'GOOD2:L,NEV2:I')
-c$$$
-c$$$ 411 format('NDBLT:I::[0,',I5,']')
-c$$$ write(block1,411) ndblt_max_nt
-c$$$ call HBNAME(ntp_level2+1,'HOUGH YZ',ndblt_nt,
-c$$$ $ block1//'
-c$$$ $ ,ALFAYZ1(NDBLT):R
-c$$$ $ ,ALFAYZ2(NDBLT):R
-c$$$ $ ,DB_CLOUD(NDBLT):I
-c$$$ $ ')
-c$$$
-c$$$ 412 format('NTRPT:I::[0,',I5,']')
-c$$$ write(block2,412) ntrpt_max_nt
-c$$$ call HBNAME(ntp_level2+1,'HOUGH XZ',NTRPT_nt,
-c$$$ $ block2//'
-c$$$ $ ,ALFAXZ1(NTRPT):R
-c$$$ $ ,ALFAXZ2(NTRPT):R
-c$$$ $ ,ALFAXZ3(NTRPT):R
-c$$$ $ ,TR_CLOUD(NTRPT):I
-c$$$ $ ')
-c$$$
-c$$$
-c$$$ 413 format('NCLOUDS_YZ:I::[0,',I4,']')
-c$$$c$$$ 414 format('DB_CLOUD(',I4,'):I')
-c$$$ write(block3,413) ncloyz_max
-c$$$c$$$ write(block4,414) ndblt_max_nt
-c$$$ call HBNAME(ntp_level2+1,'CLOUD YZ',NCLOUDS_YZ,
-c$$$ $ block3//'
-c$$$ $ ,ALFAYZ1_AV(NCLOUDS_YZ):R
-c$$$ $ ,ALFAYZ2_AV(NCLOUDS_YZ):R
-c$$$ $ ,PTCLOUD_YZ(NCLOUDS_YZ):I'
-c$$$c$$$ $ ,'//block4
-c$$$ $ )
-c$$$
-c$$$ 415 format('NCLOUDS_XZ:I::[0,',I4,']')
-c$$$c$$$ 416 format('TR_CLOUD(',I5,'):I')
-c$$$ write(block5,415) ncloxz_max
-c$$$c$$$ write(block6,416) ntrpt_max_nt
-c$$$ call HBNAME(ntp_level2+1,'CLOUD XZ',NCLOUDS_XZ,
-c$$$ $ block5//'
-c$$$ $ ,ALFAXZ1_AV(NCLOUDS_XZ):R
-c$$$ $ ,ALFAXZ2_AV(NCLOUDS_XZ):R
-c$$$ $ ,ALFAXZ3_AV(NCLOUDS_XZ):R
-c$$$ $ ,PTCLOUD_XZ(NCLOUDS_XZ):I'
-c$$$c$$$ $ ,'//block6
-c$$$ $ )
-c$$$
-c$$$
-c$$$ return
-c$$$ end
-***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...***
-*
-*
-*
-*
-*
-*
-*
-*
-*
-***...***...***...***...***...***...***...***...***...***...***...***...***...***...***...***
-c$$$ subroutine book_level2
-c$$$c*****************************************************
-c$$$cccccc 11/9/2005 modified by david fedele
-c$$$cccccc 07/10/2005 modified by elena vannuccini --> (2)
-c$$$c*****************************************************
-c$$$
-c$$$ include 'commontracker.f'
-c$$$ include 'common_momanhough.f'
-c$$$ include 'level2.f'
-c$$$
-c$$$ character*35 block1,block2
-c$$$
-c$$$c print*,'__________ booking LEVEL2 n-tuple __________'
-c$$$
-c$$$* LEVEL1 ntuple:
-c$$$ call HBNT(ntp_level2,'LEVEL2',' ')
-c$$$
-c$$$c*****************************************************
-c$$$cccccc 11/9/2005 modified by david fedele
-c$$$c call HBNAME(ntp_level2,'EVENT',good2,'GOOD2:L,NEV2:I')
-c$$$cccccc 06/10/2005 modified by elena vannuccini
-c$$$c call HBNAME(ntp_level2,'GENERAL',good2,'GOOD2:L,NEV2:I
-c$$$c $ ,WHIC_CALIB:I,SWCODE:I')
-c$$$ call HBNAME(ntp_level2,'GENERAL',good2,'GOOD2:L,NEV2:I
-c$$$ $ ,WHICH_CALIB:I,SWCODE:I,CRC(12):L')
-c$$$c*********************************************************
-c$$$
-c$$$
-c$$$c$$$# ifndef TEST2003
-c$$$c$$$
-c$$$c$$$ call HBNAME(ntp_level2,'CPU',pkt_type
-c$$$c$$$ $ ,'PKT_TYPE:I::[0,50]
-c$$$c$$$ $ ,PKT_NUM:I
-c$$$c$$$ $ ,OBT:I'//
-c$$$c$$$c********************************************************
-c$$$c$$$cccccc 11/9/2005 modified by david fedele
-c$$$c$$$c $ ,WHICH_CALIB:I::[0,50]')
-c$$$c$$$ $ ',CPU_CRC:L')
-c$$$c$$$c********************************************************
-c$$$c$$$
-c$$$c$$$# endif
-c$$$
-c$$$ 417 format('NTRK:I::[0,',I4,']')
-c$$$ 418 format(',IMAGE(NTRK):I::[0,',I4,']')
-c$$$ write(block1,417)NTRKMAX
-c$$$ write(block2,418)NTRKMAX
-c$$$ call HBNAME(ntp_level2,'TRACKS',NTRK,
-c$$$ $ block1//
-c$$$ $ block2//'
-c$$$ $ ,XM(6,NTRK):R
-c$$$ $ ,YM(6,NTRK):R
-c$$$ $ ,ZM(6,NTRK):R
-c$$$ $ ,RESX(6,NTRK):R
-c$$$ $ ,RESY(6,NTRK):R
-c$$$ $ ,AL(5,NTRK):R
-c$$$ $ ,COVAL(5,5,NTRK):R
-c$$$ $ ,CHI2(NTRK):R
-c$$$ $ ,XGOOD(6,NTRK):I::[0,1]
-c$$$ $ ,YGOOD(6,NTRK):I::[0,1]
-c$$$ $ ,XV(6,NTRK):R
-c$$$ $ ,YV(6,NTRK):R
-c$$$ $ ,ZV(6,NTRK):R
-c$$$ $ ,AXV(6,NTRK):R
-c$$$ $ ,AYV(6,NTRK):R'//
-c$$$c*****************************************************
-c$$$cccccc 11/9/2005 modified by david fedele
-c$$$c $ ,DEDXP(6,NTRK):R'//
-c$$$c $ ')
-c$$$ $ ',DEDX_X(6,NTRK):R
-c$$$ $ ,DEDX_Y(6,NTRK):R'//
-c$$$c****************************************************
-c$$$cccccc 06/10/2005 modified by elena vannuccini
-c$$$c $ ,CRC(12):L
-c$$$ $ ',BdL(NTRK):R'
-c$$$ $ )
-c$$$c****************************************************
-c$$$
-c$$$
-c$$$ call HBNAME(ntp_level2,'SINGLETX',nclsx,
-c$$$c*****************************************************
-c$$$cccccc 11/9/2005 modified by david fedele
-c$$$c $ 'NCLSX(6):I,NCLSY(6):I')
-c$$$ $ 'NCLSX:I::[0,500],PLANEX(NCLSX):I
-c$$$ $ ,XS(2,NCLSX):R,SGNLXS(NCLSX):R') !(2)
-c$$$c $ ,XS(NCLSX):R,SGNLXS(NCLSX):R') !(2)
-c$$$ call HBNAME(ntp_level2,'SINGLETY',nclsy,
-c$$$ $ 'NCLSY:I::[0,500],PLANEY(NCLSY):I
-c$$$ $ ,YS(2,NCLSY):R,SGNLYS(NCLSY):R') !(2)
-c$$$c $ ,YS(NCLSY):R,SGNLYS(NCLSY):R') !(2)
-c$$$ return
-c$$$ end
-
-c$$$ subroutine fill_level2_clouds
-c$$$c*****************************************************
-c$$$c 29/11/2005 created by elena vannuccini
-c$$$c*****************************************************
-c$$$
-c$$$* -------------------------------------------------------
-c$$$* This routine fills the variables related to the hough
-c$$$* transform, for the debig n-tuple
-c$$$* -------------------------------------------------------
-c$$$
-c$$$ include 'commontracker.f'
-c$$$ include 'common_momanhough.f'
-c$$$ include 'common_level2debug.f'
-c$$$ include 'level2.f'
-c$$$
-c$$$ good2_nt=.true.!good2
-c$$$c nev2_nt=nev2
-c$$$
-c$$$ if(.false.
-c$$$ $ .or.ntrpt.gt.ntrpt_max_nt
-c$$$ $ .or.ndblt.gt.ndblt_max_nt
-c$$$ $ .or.NCLOUDS_XZ.gt.ncloxz_max
-c$$$ $ .or.NCLOUDS_yZ.gt.ncloyz_max
-c$$$ $ )then
-c$$$ good2_nt=.false.
-c$$$ ntrpt_nt=0
-c$$$ ndblt_nt=0
-c$$$ NCLOUDS_XZ_nt=0
-c$$$ NCLOUDS_YZ_nt=0
-c$$$ else
-c$$$ ndblt_nt=ndblt
-c$$$ ntrpt_nt=ntrpt
-c$$$ if(ndblt.ne.0)then
-c$$$ do id=1,ndblt
-c$$$ alfayz1_nt(id)=alfayz1(id) !Y0
-c$$$ alfayz2_nt(id)=alfayz2(id) !tg theta-yz
-c$$$c db_cloud_nt(id)=db_cloud(id)
-c$$$ enddo
-c$$$ endif
-c$$$ if(ndblt.ne.0)then
-c$$$ do it=1,ntrpt
-c$$$ alfaxz1_nt(it)=alfaxz1(it) !X0
-c$$$ alfaxz2_nt(it)=alfaxz2(it) !tg theta-xz
-c$$$ alfaxz3_nt(it)=alfaxz3(it) !1/r
-c$$$c tr_cloud_nt(it)=tr_cloud(it)
-c$$$ enddo
-c$$$ endif
-c$$$ nclouds_yz_nt=nclouds_yz
-c$$$ nclouds_xz_nt=nclouds_xz
-c$$$ if(nclouds_yz.ne.0)then
-c$$$ nnn=0
-c$$$ do iyz=1,nclouds_yz
-c$$$ ptcloud_yz_nt(iyz)=ptcloud_yz(iyz)
-c$$$ alfayz1_av_nt(iyz)=alfayz1_av(iyz)
-c$$$ alfayz2_av_nt(iyz)=alfayz2_av(iyz)
-c$$$ nnn=nnn+ptcloud_yz(iyz)
-c$$$ enddo
-c$$$ do ipt=1,nnn
-c$$$ db_cloud_nt(ipt)=db_cloud(ipt)
-c$$$ enddo
-c$$$c print*,'#### ntupla #### ptcloud_yz '
-c$$$c $ ,(ptcloud_yz(i),i=1,nclouds_yz)
-c$$$c print*,'#### ntupla #### db_cloud '
-c$$$c $ ,(db_cloud(i),i=1,nnn)
-c$$$ endif
-c$$$ if(nclouds_xz.ne.0)then
-c$$$ nnn=0
-c$$$ do ixz=1,nclouds_xz
-c$$$ ptcloud_xz_nt(ixz)=ptcloud_xz(ixz)
-c$$$ alfaxz1_av_nt(ixz)=alfaxz1_av(ixz)
-c$$$ alfaxz2_av_nt(ixz)=alfaxz2_av(ixz)
-c$$$ alfaxz3_av_nt(ixz)=alfaxz3_av(ixz)
-c$$$ nnn=nnn+ptcloud_xz(ixz)
-c$$$ enddo
-c$$$ do ipt=1,nnn
-c$$$ tr_cloud_nt(ipt)=tr_cloud(ipt)
-c$$$ enddo
-c$$$c print*,'#### ntupla #### ptcloud_xz '
-c$$$c $ ,(ptcloud_xz(i),i=1,nclouds_xz)
-c$$$c print*,'#### ntupla #### tr_cloud '
-c$$$c $ ,(tr_cloud(i),i=1,nnn)
-c$$$ endif
-c$$$ endif
-c$$$ end
***************************************************
@@ -1726,8 +1595,12 @@
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
- ncp_tot=0
do ip=1,nplanes
do ico=1,ncouplemax
clx(ip,ico)=0
@@ -1740,8 +1613,8 @@
ncls(ip)=0
enddo
do icl=1,nclstrmax_level2
- cl_single(icl) = 1
- cl_good(icl) = 0
+ cl_single(icl) = 1 !all are single
+ cl_good(icl) = 0 !all are bad
enddo
do iv=1,nviews
ncl_view(iv) = 0
@@ -1754,10 +1627,12 @@
enddo
* mask views with too many clusters
do iv=1,nviews
- if( ncl_view(iv).gt. nclustermax)then
- mask_view(iv) = 1
- if(VERBOSE)print*,' * WARNING * cl_to_couple: n.clusters > '
- $ ,nclustermax,' on view ', iv,' --> masked!'
+ 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
@@ -1767,6 +1642,8 @@
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)
* ----------------------------------------------------
@@ -1774,7 +1651,7 @@
* ----------------------------------------------------
* cut on charge (X VIEW)
* ----------------------------------------------------
- if(dedx(icx).lt.dedx_x_min)then
+ if(sgnl(icx)/mip(VIEW(icx),LADDER(icx)).lt.dedx_x_min)then
cl_single(icx)=0
goto 10
endif
@@ -1815,6 +1692,8 @@
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)
@@ -1824,7 +1703,7 @@
* ----------------------------------------------------
* cut on charge (Y VIEW)
* ----------------------------------------------------
- if(dedx(icy).lt.dedx_y_min)then
+ if(sgnl(icy)/mip(VIEW(icy),LADDER(icy)) .lt.dedx_y_min)then
cl_single(icy)=0
goto 20
endif
@@ -1870,21 +1749,21 @@
* charge correlation
* (modified to be applied only below saturation... obviously)
- if( .not.(dedx(icy).gt.chsaty.and.dedx(icx).gt.chsatx)
+ if( .not.(sgnl(icy).gt.chsaty.and.sgnl(icx).gt.chsatx)
$ .and.
- $ .not.(dedx(icy).lt.chmipy.and.dedx(icx).lt.chmipx)
+ $ .not.(sgnl(icy).lt.chmipy.and.sgnl(icx).lt.chmipx)
$ .and.
$ (badclx.eq.1.and.badcly.eq.1)
$ .and.
$ .true.)then
- ddd=(dedx(icy)
- $ -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))
+ 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)*dedx(icy)+dedx(icx)
+ 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.)
@@ -1893,50 +1772,36 @@
goto 20 !charge not consistent
endif
endif
-
-* ------------------> COUPLE <------------------
-* check to do not overflow vector dimentions
-c$$$ if(ncp_plane(nplx).gt.ncouplemax)then
-c$$$ if(DEBUG)print*,
-c$$$ $ ' ** warning ** number of identified'//
-c$$$ $ ' couples on plane ',nplx,
-c$$$ $ ' exceeds vector dimention'//
-c$$$ $ ' ( ',ncouplemax,' )'
-c$$$c good2=.false.
-c$$$c goto 880 !fill ntp and go to next event
-c$$$ iflag=1
-c$$$ return
-c$$$ endif
-
- if(ncp_plane(nplx).eq.ncouplemax)then
- if(verbose)print*,
+
+ if(ncp_plane(nplx).gt.ncouplemax)then
+ if(verbose.eq.1)print*,
$ '** warning ** number of identified '//
$ 'couples on plane ',nplx,
$ 'exceeds vector dimention '
- $ ,'( ',ncouplemax,' )'
-c good2=.false.
-c goto 880 !fill ntp and go to next event
- mask_view(nviewx(nplx)) = 2
- mask_view(nviewy(nply)) = 2
-c iflag=1
-c return
+ $ ,'( ',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
* ----------------------------------------------
+ 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))
@@ -1944,28 +1809,86 @@
cls(ip,ncls(ip))=icl
endif
enddo
+
+c 80 continue
+ continue
- if(DEBUG)then
+ if(DEBUG.EQ.1)then
print*,'clusters ',nclstr1
print*,'good ',(cl_good(i),i=1,nclstr1)
- print*,'singles ',(cl_single(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
- do ip=1,6
+ ncp_tot=0
+ do ip=1,NPLANES
ncp_tot = ncp_tot + ncp_plane(ip)
enddo
-c if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)
-
-c$$$ if(ncp_tot.gt.ncp_max)then
-c$$$ if(verbose)print*,
-c$$$ $ '** warning ** number of identified '//
-c$$$ $ 'couples exceeds upper limit for Hough tr. '
-c$$$ $ ,'( ',ncp_max,' )'
-c$$$ iflag=1
-c$$$ return
-c$$$ endif
+ if(DEBUG.EQ.1)
+ $ print*,'n.couple tot: ',ncp_tot
return
end
@@ -1978,240 +1901,15 @@
* *
* *
**************************************************
-c$$$ subroutine cl_to_couples_nocharge(iflag)
-c$$$
-c$$$ include 'commontracker.f'
-c$$$ include 'level1.f'
-c$$$ include 'common_momanhough.f'
-c$$$c include 'momanhough_init.f'
-c$$$ include 'calib.f'
-c$$$c include 'level1.f'
-c$$$
-c$$$
-c$$$* output flag
-c$$$* --------------
-c$$$* 0 = good event
-c$$$* 1 = bad event
-c$$$* --------------
-c$$$ integer iflag
-c$$$
-c$$$ integer badseed,badcl
-c$$$
-c$$$* init variables
-c$$$ ncp_tot=0
-c$$$ do ip=1,nplanes
-c$$$ do ico=1,ncouplemax
-c$$$ clx(ip,ico)=0
-c$$$ cly(ip,ico)=0
-c$$$ enddo
-c$$$ ncp_plane(ip)=0
-c$$$ do icl=1,nclstrmax_level2
-c$$$ cls(ip,icl)=1
-c$$$ enddo
-c$$$ ncls(ip)=0
-c$$$ enddo
-c$$$ do icl=1,nclstrmax_level2
-c$$$ cl_single(icl)=1
-c$$$ cl_good(icl)=0
-c$$$ enddo
-c$$$
-c$$$* start association
-c$$$ ncouples=0
-c$$$ do icx=1,nclstr1 !loop on cluster (X)
-c$$$ if(mod(VIEW(icx),2).eq.1)goto 10
-c$$$
-c$$$* ----------------------------------------------------
-c$$$* cut on charge (X VIEW)
-c$$$ if(dedx(icx).lt.dedx_x_min)then
-c$$$ cl_single(icx)=0
-c$$$ goto 10
-c$$$ endif
-c$$$* cut BAD (X VIEW)
-c$$$ badseed=BAD(VIEW(icx),nvk(MAXS(icx)),nst(MAXS(icx)))
-c$$$ ifirst=INDSTART(icx)
-c$$$ if(icx.ne.nclstr1) then
-c$$$ ilast=INDSTART(icx+1)-1
-c$$$ else
-c$$$ ilast=TOTCLLENGTH
-c$$$ endif
-c$$$ badcl=badseed
-c$$$ do igood=-ngoodstr,ngoodstr
-c$$$ ibad=1
-c$$$ if((INDMAX(icx)+igood).gt.ifirst.and.
-c$$$ $ (INDMAX(icx)+igood).lt.ilast.and.
-c$$$ $ .true.)then
-c$$$ ibad=BAD(VIEW(icx),
-c$$$ $ nvk(MAXS(icx)+igood),
-c$$$ $ nst(MAXS(icx)+igood))
-c$$$ endif
-c$$$ badcl=badcl*ibad
-c$$$ enddo
-c$$$ if(badcl.eq.0)then !<<<<<<<<<<<<<< BAD cut
-c$$$ cl_single(icx)=0 !<<<<<<<<<<<<<< BAD cut
-c$$$ goto 10 !<<<<<<<<<<<<<< BAD cut
-c$$$ endif !<<<<<<<<<<<<<< BAD cut
-c$$$* ----------------------------------------------------
-c$$$
-c$$$ cl_good(icx)=1
-c$$$ nplx=npl(VIEW(icx))
-c$$$ nldx=nld(MAXS(icx),VIEW(icx))
-c$$$
-c$$$ do icy=1,nclstr1 !loop on cluster (Y)
-c$$$ if(mod(VIEW(icy),2).eq.0)goto 20
-c$$$
-c$$$* ----------------------------------------------------
-c$$$* cut on charge (Y VIEW)
-c$$$ if(dedx(icy).lt.dedx_y_min)then
-c$$$ cl_single(icy)=0
-c$$$ goto 20
-c$$$ endif
-c$$$* cut BAD (Y VIEW)
-c$$$ badseed=BAD(VIEW(icy),nvk(MAXS(icy)),nst(MAXS(icy)))
-c$$$ ifirst=INDSTART(icy)
-c$$$ if(icy.ne.nclstr1) then
-c$$$ ilast=INDSTART(icy+1)-1
-c$$$ else
-c$$$ ilast=TOTCLLENGTH
-c$$$ endif
-c$$$ badcl=badseed
-c$$$ do igood=-ngoodstr,ngoodstr
-c$$$ ibad=1
-c$$$ if((INDMAX(icy)+igood).gt.ifirst.and.
-c$$$ $ (INDMAX(icy)+igood).lt.ilast.and.
-c$$$ $ .true.)
-c$$$ $ ibad=BAD(VIEW(icy),
-c$$$ $ nvk(MAXS(icy)+igood),
-c$$$ $ nst(MAXS(icy)+igood))
-c$$$ badcl=badcl*ibad
-c$$$ enddo
-c$$$ if(badcl.eq.0)then !<<<<<<<<<<<<<< BAD cut
-c$$$ cl_single(icy)=0 !<<<<<<<<<<<<<< BAD cut
-c$$$ goto 20 !<<<<<<<<<<<<<< BAD cut
-c$$$ endif !<<<<<<<<<<<<<< BAD cut
-c$$$* ----------------------------------------------------
-c$$$
-c$$$
-c$$$ cl_good(icy)=1
-c$$$ nply=npl(VIEW(icy))
-c$$$ nldy=nld(MAXS(icy),VIEW(icy))
-c$$$
-c$$$* ----------------------------------------------
-c$$$* CONDITION TO FORM A COUPLE
-c$$$* ----------------------------------------------
-c$$$* geometrical consistency (same plane and ladder)
-c$$$ if(nply.eq.nplx.and.nldy.eq.nldx)then
-c$$$* charge correlation
-c$$$* ===========================================================
-c$$$* this version of the subroutine is used for the calibration
-c$$$* thus charge-correlation selection is obviously removed
-c$$$* ===========================================================
-c$$$c$$$ ddd=(dedx(icy)
-c$$$c$$$ $ -kch(nplx,nldx)*dedx(icx)-cch(nplx,nldx))
-c$$$c$$$ ddd=ddd/sqrt(kch(nplx,nldx)**2+1)
-c$$$c$$$ cut=chcut*sch(nplx,nldx)
-c$$$c$$$ if(abs(ddd).gt.cut)goto 20 !charge not consistent
-c$$$* ===========================================================
-c$$$
-c$$$
-c$$$* ------------------> COUPLE <------------------
-c$$$* check to do not overflow vector dimentions
-c$$$c$$$ if(ncp_plane(nplx).gt.ncouplemax)then
-c$$$c$$$ if(DEBUG)print*,
-c$$$c$$$ $ ' ** warning ** number of identified'//
-c$$$c$$$ $ ' couples on plane ',nplx,
-c$$$c$$$ $ ' exceeds vector dimention'//
-c$$$c$$$ $ ' ( ',ncouplemax,' )'
-c$$$c$$$c good2=.false.
-c$$$c$$$c goto 880 !fill ntp and go to next event
-c$$$c$$$ iflag=1
-c$$$c$$$ return
-c$$$c$$$ endif
-c$$$
-c$$$ if(ncp_plane(nplx).eq.ncouplemax)then
-c$$$ if(verbose)print*,
-c$$$ $ '** warning ** number of identified '//
-c$$$ $ 'couples on plane ',nplx,
-c$$$ $ 'exceeds vector dimention '
-c$$$ $ ,'( ',ncouplemax,' )'
-c$$$c good2=.false.
-c$$$c goto 880 !fill ntp and go to next event
-c$$$ iflag=1
-c$$$ return
-c$$$ endif
-c$$$
-c$$$ ncp_plane(nplx) = ncp_plane(nplx) + 1
-c$$$ clx(nplx,ncp_plane(nplx))=icx
-c$$$ cly(nply,ncp_plane(nplx))=icy
-c$$$ cl_single(icx)=0
-c$$$ cl_single(icy)=0
-c$$$ endif
-c$$$* ----------------------------------------------
-c$$$
-c$$$ 20 continue
-c$$$ enddo !end loop on clusters(Y)
-c$$$
-c$$$ 10 continue
-c$$$ enddo !end loop on clusters(X)
-c$$$
-c$$$
-c$$$ do icl=1,nclstr1
-c$$$ if(cl_single(icl).eq.1)then
-c$$$ ip=npl(VIEW(icl))
-c$$$ ncls(ip)=ncls(ip)+1
-c$$$ cls(ip,ncls(ip))=icl
-c$$$ endif
-c$$$ enddo
-c$$$
-c$$$
-c$$$ if(DEBUG)then
-c$$$ print*,'clusters ',nclstr1
-c$$$ print*,'good ',(cl_good(i),i=1,nclstr1)
-c$$$ print*,'singles ',(cl_single(i),i=1,nclstr1)
-c$$$ print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)
-c$$$ endif
-c$$$
-c$$$ do ip=1,6
-c$$$ ncp_tot=ncp_tot+ncp_plane(ip)
-c$$$ enddo
-c$$$c if(ncp_tot.gt.ncp_max)goto 100!next event (TEMPORANEO!!!)
-c$$$
-c$$$ if(ncp_tot.gt.ncp_max)then
-c$$$ if(verbose)print*,
-c$$$ $ '** warning ** number of identified '//
-c$$$ $ 'couples exceeds upper limit for Hough tr. '
-c$$$ $ ,'( ',ncp_max,' )'
-c$$$c good2=.false.
-c$$$c goto 880 !fill ntp and go to next event
-c$$$ iflag=1
-c$$$ return
-c$$$ endif
-c$$$
-c$$$ return
-c$$$ end
-c$$$
-
-***************************************************
-* *
-* *
-* *
-* *
-* *
-* *
-**************************************************
subroutine cp_to_doubtrip(iflag)
-c*****************************************************
-c 02/02/2006 modified by Elena Vannuccini --> (1)
-c*****************************************************
include 'commontracker.f'
include 'level1.f'
include 'common_momanhough.f'
-c include 'momanhough_init.f'
include 'common_xyzPAM.f'
include 'common_mini_2.f'
include 'calib.f'
-c include 'level1.f'
* output flag
@@ -2244,54 +1942,99 @@
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
- do is1=1,2 !loop on sensors - 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)
- call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)
- xm1=xPAM
- ym1=yPAM
- zm1=zPAM
-c print*,'***',is1,xm1,ym1,zm1
+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
- do is2=1,2 !loop on sensors -ndblt 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.) !(1)
- xm2=xPAM
- ym2=yPAM
- zm2=zPAM
-
+ $ (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)print*,
+ 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
+c good2=.false.
+c goto 880 !fill ntp and go to next event
do iv=1,12
- mask_view(iv) = 3
+c mask_view(iv) = 3
+ mask_view(iv) = mask_view(iv)+ 2**2
enddo
iflag=1
return
endif
+
+
+ccc print*,'<doublet> ',icp1,icp2
+
ndblt = ndblt + 1
* store doublet info
cpyz1(ndblt)=id_cp(ip1,icp1,is1)
@@ -2299,119 +2042,205 @@
* tg(th_yz)
alfayz2(ndblt)=(ym1-ym2)/(zm1-zm2)
* y0 (cm)
- alfayz1(ndblt)=alfayz2(ndblt)*(zini-zm1)+ym1
-
+ 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
-c$$$ if(iev.eq.33)then
-c$$$ print*,'********* ',ndblt,' -- ',icp1,icp2,is1,is2
-c$$$ $ ,' || ',icx1,icy1,icx2,icy2
-c$$$ $ ,' || ',xm1,ym1,xm2,ym2
-c$$$ $ ,' || ',alfayz2(ndblt),alfayz1(ndblt)
-c$$$ endif
-c$$$
+
+ 111 continue
* - - - - - - - - - - - - - - - - - - - - - - - - - - - -
* track parameters on Y VIEW - end
* - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- if(ip2.eq.nplanes)goto 30 !no possible combination with 3 couples
+ 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.) !(1)
- xm3=xPAM
- ym3=yPAM
- zm3=zPAM
+ $ (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)
-c print*,xc,zc,radius
+ $ ,xc,zc,radius,iflag_t)
* the circle must intersect the reference plane
- if(
-c $ (xc.le.-1.*xclimit.or.
-c $ xc.ge.xclimit).and.
- $ radius**2.ge.(ZINI-zc)**2.and.
- $ iflag.eq.0.and.
- $ .true.)then
-
+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)print*,
- $ '** warning ** number of identified '//
- $ 'triplets exceeds vector dimention '
- $ ,'( ',ntrpt_max,' )'
-c good2=.false.
-c goto 880 !fill ntp and go to next event
+ 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
- mask_view(iv) = 4
+c mask_view(iv) = 4
+ mask_view(iv) =
+ $ mask_view(iv)+ 2**3
enddo
iflag=1
return
endif
+
+ccc print*,'<triplet> ',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(xc.lt.0)then
+ if(radius.ne.0.and.xc.lt.0)then
*************POSITIVE DEFLECTION
- alfaxz1(ntrpt) = xc+sqrt(radius**2-(ZINI-zc)**2)
- alfaxz2(ntrpt) = (ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)
- alfaxz3(ntrpt) = 1/radius
- else
+ 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-(ZINI-zc)**2)
- alfaxz2(ntrpt) = -(ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)
- alfaxz3(ntrpt) = -1/radius
- endif
-
+ 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(
- $ abs(alfaxz2(ntrpt)).gt.alfxz2_max
- $ .or.
- $ abs(alfaxz1(ntrpt)).gt.alfxz1_max
- $ )ntrpt = ntrpt-1
-
-
-c print*,alfaxz1(ntrpt),alfaxz2(ntrpt),alfaxz3(ntrpt)
+ 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
* - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- endif
+
+ 29 continue
enddo !end loop on COPPIA 3
enddo !end loop on sensors - COPPIA 3
+ 30 continue
enddo !end loop on planes - COPPIA 3
- 30 continue
-
- 1 enddo !end loop on COPPIA 2
+
+ 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)then
+ if(DEBUG.EQ.1)then
print*,'--- doublets ',ndblt
print*,'--- triplets ',ntrpt
endif
@@ -2458,6 +2287,7 @@
integer cp_useds1(ncouplemaxtot) ! sensor 1
integer cp_useds2(ncouplemaxtot) ! sensor 2
+ if(DEBUG.EQ.1)print*,'doub_to_YZcloud:'
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
* classification of DOUBLETS
@@ -2476,9 +2306,6 @@
do idb1=1,ndblt !loop (1) on DOUBLETS
if(db_used(idb1).eq.1)goto 2228 !db already included in a cloud
-c print*,'--------------'
-c print*,'** ',idb1,' **'
-
do icp=1,ncp_tot
cp_useds1(icp)=0 !init
cp_useds2(icp)=0 !init
@@ -2514,19 +2341,11 @@
* doublet distance in parameter space
distance=
$ ((alfayz1(idbref)-alfayz1(idb2))/Dalfayz1)**2
- $ +((alfayz2(idbref)-alfayz2(idb2))/Dalfayz2)**2
+ $ +((alfayz2(idbref)-alfayz2(idb2))/Dalfayz2)**2
distance = sqrt(distance)
-c$$$ if(iev.eq.33)then
-c$$$ if(distance.lt.100)
-c$$$ $ print*,'********* ',idb1,idbref,idb2,distance
-c$$$ if(distance.lt.100)
-c$$$ $ print*,'********* ',alfayz1(idbref),alfayz1(idb2)
-c$$$ $ ,alfayz2(idbref),alfayz2(idb2)
-c$$$ endif
if(distance.lt.cutdistyz)then
-c print*,idb1,idb2,distance,' cloud ',nclouds_yz
if(cpyz1(idb2).gt.0)cp_useds2(cpyz1(idb2))=1
if(cpyz1(idb2).lt.0)cp_useds1(-cpyz1(idb2))=1
if(cpyz2(idb2).gt.0)cp_useds2(cpyz2(idb2))=1
@@ -2542,13 +2361,13 @@
temp1 = temp1 + alfayz1(idb2)
temp2 = temp2 + alfayz2(idb2)
-c print*,'* idbref,idb2 ',idbref,idb2
endif
1118 continue
enddo !end loop (2) on DOUBLETS
- 1188 continue
+c 1188 continue
+ continue
enddo !end loop on... bo?
nptloop=npv
@@ -2565,7 +2384,9 @@
enddo
ncpused=0
do icp=1,ncp_tot
- if(cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0)then
+ 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
@@ -2575,23 +2396,22 @@
do ip=1,nplanes
nplused=nplused+ hit_plane(ip)
enddo
-c print*,'>>>> ',ncpused,npt,nplused
-c if(ncpused.lt.ncpyz_min)goto 2228 !next doublet
- if(npt.lt.nptyz_min)goto 2228 !next doublet
+
if(nplused.lt.nplyz_min)goto 2228 !next doublet
* ~~~~~~~~~~~~~~~~~
* >>> NEW CLOUD <<<
if(nclouds_yz.ge.ncloyz_max)then
- if(verbose)print*,
+ 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
- mask_view(iv) = 5
+c mask_view(iv) = 5
+ mask_view(iv) = mask_view(iv) + 2**4
enddo
iflag=1
return
@@ -2608,21 +2428,17 @@
c ptcloud_yz_nt(nclouds_yz)=npt
do ipt=1,npt
db_cloud(npt_tot+ipt) = db_all(ipt)
-c print*,'>> ',ipt,db_all(ipt)
enddo
npt_tot=npt_tot+npt
- if(DEBUG)then
- print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'
+ 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*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)
-c$$$ print*,'nt-uple: ptcloud_yz(',nclouds_yz,') = '
-c$$$ $ ,ptcloud_yz(nclouds_yz)
-c$$$ print*,'nt-uple: db_cloud(...) = '
-c$$$ $ ,(db_cloud(iii),iii=npt_tot-npt+1,npt_tot)
+ 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 <<<
* ~~~~~~~~~~~~~~~~~
@@ -2636,10 +2452,8 @@
goto 90
endif
- if(DEBUG)then
- print*,'---------------------- '
+ if(DEBUG.EQ.1)then
print*,'Y-Z total clouds ',nclouds_yz
- print*,' '
endif
@@ -2685,6 +2499,8 @@
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
@@ -2700,8 +2516,6 @@
91 continue
do itr1=1,ntrpt !loop (1) on TRIPLETS
if(tr_used(itr1).eq.1)goto 22288 !already included in a cloud
-c print*,'--------------'
-c print*,'** ',itr1,' **'
do icp=1,ncp_tot
cp_useds1(icp)=0
@@ -2735,15 +2549,28 @@
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
+ $ +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2
distance = sqrt(distance)
-
- if(distance.lt.cutdistxz)then
-c print*,idb1,idb2,distance,' cloud ',nclouds_yz
+
+
+* ------------------------------------------------------------------------
+* 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
@@ -2762,13 +2589,13 @@
temp1 = temp1 + alfaxz1(itr2)
temp2 = temp2 + alfaxz2(itr2)
temp3 = temp3 + alfaxz3(itr2)
-c print*,'* itrref,itr2 ',itrref,itr2,distance
endif
11188 continue
enddo !end loop (2) on TRIPLETS
-11888 continue
+c11888 continue
+ continue
enddo !end loop on... bo?
nptloop=npv
@@ -2783,13 +2610,14 @@
* 1bis)
* 2) it is not already stored
* ------------------------------------------
-c print*,'check cp_used'
do ip=1,nplanes
hit_plane(ip)=0
enddo
ncpused=0
do icp=1,ncp_tot
- if(cp_useds1(icp).ne.0.or.cp_useds2(icp).ne.0)then
+ 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
@@ -2799,21 +2627,20 @@
do ip=1,nplanes
nplused=nplused+ hit_plane(ip)
enddo
-c if(ncpused.lt.ncpxz_min)goto 22288 !next triplet
- if(npt.lt.nptxz_min)goto 22288 !next triplet
- if(nplused.lt.nplxz_min)goto 22288 !next doublet
+ if(nplused.lt.nplxz_min)goto 22288 !next triplet
* ~~~~~~~~~~~~~~~~~
* >>> NEW CLOUD <<<
if(nclouds_xz.ge.ncloxz_max)then
- if(verbose)print*,
+ 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
- mask_view(iv) = 6
+c mask_view(iv) = 6
+ mask_view(iv) = mask_view(iv) + 2**5
enddo
iflag=1
return
@@ -2832,19 +2659,16 @@
enddo
npt_tot=npt_tot+npt
- if(DEBUG)then
- print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'
- print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'
- print*,'- alfaxz1 ',alfaxz1_av(nclouds_xz)
- print*,'- alfaxz2 ',alfaxz2_av(nclouds_xz)
- print*,'- alfaxz3 ',alfaxz3_av(nclouds_xz)
- print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot)
- print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot)
+ 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)
-c$$$ print*,'nt-uple: ptcloud_xz(',nclouds_xz,') = '
-c$$$ $ ,ptcloud_xz(nclouds_xz)
-c$$$ print*,'nt-uple: tr_cloud(...) = '
-c$$$ $ ,(tr_cloud(iii),iii=npt_tot-npt+1,npt_tot)
endif
* >>> NEW CLOUD <<<
* ~~~~~~~~~~~~~~~~~
@@ -2857,10 +2681,8 @@
goto 91
endif
- if(DEBUG)then
- print*,'---------------------- '
+ if(DEBUG.EQ.1)then
print*,'X-Z total clouds ',nclouds_xz
- print*,' '
endif
@@ -2878,9 +2700,6 @@
**************************************************
subroutine clouds_to_ctrack(iflag)
-c*****************************************************
-c 02/02/2006 modified by Elena Vannuccini --> (1)
-c*****************************************************
include 'commontracker.f'
include 'level1.f'
@@ -2888,7 +2707,7 @@
include 'common_xyzPAM.f'
include 'common_mini_2.f'
include 'common_mech.f'
-c include 'momanhough_init.f'
+
* output flag
@@ -2912,16 +2731,15 @@
* -----------------------------------------------------------
* variables for track fitting
double precision AL_INI(5)
- double precision tath
* -----------------------------------------------------------
-c real fitz(nplanes) !z coordinates of the planes in cm
+ if(DEBUG.EQ.1)print*,'clouds_to_ctrack:'
ntracks=0 !counter of track candidates
- do iyz=1,nclouds_yz !loop on YZ couds
- do ixz=1,nclouds_xz !loop on XZ couds
+ do iyz=1,nclouds_yz !loop on YZ clouds
+ do ixz=1,nclouds_xz !loop on XZ clouds
* --------------------------------------------------
* check of consistency of the clouds
@@ -2930,27 +2748,72 @@
* of the two clouds
* --------------------------------------------------
do ip=1,nplanes
- hit_plane(ip)=0
- ncp_match(ip)=0
+ 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
- do icp=1,ncp_tot !loop on couples
-* get info on
- cpintersec(icp)=min(
- $ cpcloud_yz(iyz,icp),
- $ cpcloud_xz(ixz,icp))
- if(
- $ (cpcloud_yz(iyz,icp).eq.1.and.cpcloud_xz(ixz,icp).eq.2).or.
- $ (cpcloud_yz(iyz,icp).eq.2.and.cpcloud_xz(ixz,icp).eq.1).or.
- $ .false.)cpintersec(icp)=0
+ 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
- ncp_ok=ncp_ok+1
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
@@ -2977,52 +2840,33 @@
do ip=1,nplanes
nplused=nplused+ hit_plane(ip)
enddo
-
- if(nplused.lt.nplxz_min)goto 888 !next doublet
- if(ncp_ok.lt.ncpok)goto 888 !next cloud
-
- if(DEBUG)then
- print*,'Combination ',iyz,ixz
- $ ,' db ',ptcloud_yz(iyz)
- $ ,' tr ',ptcloud_xz(ixz)
- $ ,' -----> # matching couples ',ncp_ok
- endif
-c$$$ print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'
-c$$$ print*,'Configurazione cluster XZ'
-c$$$ print*,'1 -- ',(clx(1,i),i=1,ncp_plane(1))
-c$$$ print*,'2 -- ',(clx(2,i),i=1,ncp_plane(1))
-c$$$ print*,'3 -- ',(clx(3,i),i=1,ncp_plane(1))
-c$$$ print*,'4 -- ',(clx(4,i),i=1,ncp_plane(1))
-c$$$ print*,'5 -- ',(clx(5,i),i=1,ncp_plane(1))
-c$$$ print*,'6 -- ',(clx(6,i),i=1,ncp_plane(1))
-c$$$ print*,'Configurazione cluster YZ'
-c$$$ print*,'1 -- ',(cly(1,i),i=1,ncp_plane(1))
-c$$$ print*,'2 -- ',(cly(2,i),i=1,ncp_plane(1))
-c$$$ print*,'3 -- ',(cly(3,i),i=1,ncp_plane(1))
-c$$$ print*,'4 -- ',(cly(4,i),i=1,ncp_plane(1))
-c$$$ print*,'5 -- ',(cly(5,i),i=1,ncp_plane(1))
-c$$$ print*,'6 -- ',(cly(6,i),i=1,ncp_plane(1))
-c$$$ print*,'~~~~~~~~~~~~~~~~~~~~~~~~~'
-
-* -------> INITIAL GUESS <-------
- AL_INI(1) = dreal(alfaxz1_av(ixz))
- AL_INI(2) = dreal(alfayz1_av(iyz))
- AL_INI(4) = PIGR + datan(dreal(alfayz2_av(iyz))
- $ /dreal(alfaxz2_av(ixz)))
- tath = -dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))
- AL_INI(3) = tath/sqrt(1+tath**2)
- AL_INI(5) = (1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.
-
-c print*,'*******',AL_INI(5)
- if(AL_INI(5).gt.defmax)goto 888 !next cloud
-
-c print*,'alfaxz2, alfayz2 '
-c $ ,alfaxz2_av(ixz),alfayz2_av(iyz)
-
-* -------> INITIAL GUESS <-------
-c print*,'AL_INI ',(al_ini(i),i=1,5)
-
- if(DEBUG)then
+
+ 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))
@@ -3055,10 +2899,45 @@
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
+ 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)
@@ -3072,58 +2951,105 @@
* *************************
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.) !(1)
+ $ PFAdef,PFAdef,0.,0.,0.,0.)
* *************************
* -----------------------------
- 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(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(i)=AL_INI(i)
+ AL_INI(i)=AL(i)
enddo
ifail=0 !error flag in chi^2 computation
jstep=0 !number of minimization steps
iprint=0
- if(DEBUG)iprint=1
+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)then
+ if(DEBUG.EQ.1)then
print *,
$ '*** MINIMIZATION FAILURE *** '
- $ //'(mini2 in clouds_to_ctrack)'
+ $ //'(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
- chi2=-chi2
+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)print*,
+ 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
- mask_view(iv) = 7
+c mask_view(iv) = 7
+ mask_view(iv) = mask_view(iv) + 2**6
enddo
iflag=1
return
@@ -3131,14 +3057,11 @@
ntracks = ntracks + 1
-c$$$ ndof=0
- do ip=1,nplanes
-c$$$ ndof=ndof
-c$$$ $ +int(xgood(ip))
-c$$$ $ +int(ygood(ip))
+ 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))
+ 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))
@@ -3151,23 +3074,39 @@
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
- CLS_STORE(nplanes-ip+1,ntracks)=0
+ 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
-c$$$ * Number of Degree Of Freedom
-c$$$ ndof=ndof-5
-c$$$ * reduced chi^2
-c$$$ rchi2=chi2/dble(ndof)
- RCHI2_STORE(ntracks)=chi2
+ RCHI2_STORE(ntracks)=REAL(chi2)
* --------------------------------
* STORE candidate TRACK INFO - end
@@ -3187,17 +3126,23 @@
if(ntracks.eq.0)then
iflag=1
- return
+cc return
endif
- if(DEBUG)then
- print*,'****** TRACK CANDIDATES ***********'
- print*,'# R. chi2 RIG'
- do i=1,ntracks
- print*,i,' --- ',rchi2_store(i),' --- '
- $ ,1./abs(AL_STORE(5,i))
- enddo
- print*,'***********************************'
+ 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
@@ -3216,11 +3161,6 @@
subroutine refine_track(ibest)
-c******************************************************
-cccccc 06/10/2005 modified by elena vannuccini ---> (1)
-cccccc 31/01/2006 modified by elena vannuccini ---> (2)
-cccccc 12/08/2006 modified by elena vannucicni ---> (3)
-c******************************************************
include 'commontracker.f'
include 'level1.f'
@@ -3228,15 +3168,26 @@
include 'common_xyzPAM.f'
include 'common_mini_2.f'
include 'common_mech.f'
-c include 'momanhough_init.f'
-c include 'level1.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
@@ -3245,12 +3196,60 @@
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
* -------------------------------------------------
- if(XGOOD_STORE(nplanes-ip+1,ibest).eq.1..and.
+* |||||||||||||||||||||||||||||||||||||||||||||||||
+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)
@@ -3263,47 +3262,87 @@
$ ,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,
-c $ 'ETA2','ETA2',
$ PFA,PFA,
$ AXV_STORE(nplanes-ip+1,ibest),
- $ AYV_STORE(nplanes-ip+1,ibest))
-c$$$ call xyz_PAM(icx,icy,is,
-c$$$ $ 'COG2','COG2',
-c$$$ $ 0.,
-c$$$ $ 0.)
- xm(nplanes-ip+1) = xPAM
- ym(nplanes-ip+1) = yPAM
- zm(nplanes-ip+1) = zPAM
- xgood(nplanes-ip+1) = 1
- ygood(nplanes-ip+1) = 1
- resx(nplanes-ip+1) = resxPAM
- resy(nplanes-ip+1) = resyPAM
-
-c dedxtrk(nplanes-ip+1) = (dedx(icx)+dedx(icy))/2. !(1)
- dedxtrk_x(nplanes-ip+1)=dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
- dedxtrk_y(nplanes-ip+1)=dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
+ $ 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
- xP=XV_STORE(nplanes-ip+1,ibest)
- yP=YV_STORE(nplanes-ip+1,ibest)
- zP=ZV_STORE(nplanes-ip+1,ibest)
call whichsensor(ip,xP,yP,nldt,ist)
* if the track hit the plane in a dead area, go to the next plane
if(nldt.eq.0.or.ist.eq.0)goto 133
+
+ SENSOR_STORE(nplanes-ip+1,IBEST)=ist
+ LADDER_STORE(nplanes-ip+1,IBEST)=nldt
* --------------------------------------------------------------
- if(DEBUG)then
+ if(DEBUG.EQ.1)then
print*,
$ '------ Plane ',ip,' intersected on LADDER ',nldt
$ ,' SENSOR ',ist
@@ -3314,8 +3353,7 @@
* ===========================================
* STEP 1 >>>>>>> try to include a new couple
* ===========================================
-c if(DEBUG)print*,'>>>> try to include a new couple'
- distmin=1000000.
+ distmin=100000000.
xmm = 0.
ymm = 0.
zmm = 0.
@@ -3328,23 +3366,28 @@
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 !(3)
- $ cl_used(icy).ne.0.or. !or the Y cluster is already used !(3)
+ $ 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,
-c $ 'ETA2','ETA2',
$ AXV_STORE(nplanes-ip+1,ibest),
- $ AYV_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)print*,'( couple ',id
- $ ,' ) normalized distance ',distance
+ if(DEBUG.EQ.1)
+ $ print*,'( couple ',id
+ $ ,' ) distance ',distance
if(distance.lt.distmin)then
xmm = xPAM
ymm = yPAM
@@ -3353,35 +3396,34 @@
rymm = resyPAM
distmin = distance
idm = id
-c dedxmm = (dedx(icx)+dedx(icy))/2. !(1)
- dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
- dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
+ 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.clinc)then
+ 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 !<<<
-c dedxtrk(nplanes-ip+1) = dedxmm !<<< !(1)
- dedxtrk_x(nplanes-ip+1) = dedxmmx !(1)
- dedxtrk_y(nplanes-ip+1) = dedxmmy !(1)
+ 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)print*,'%%%% included couple ',idm
- $ ,' (norm.dist.= ',distmin,', cut ',clinc,' )'
+ 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
* ================================================
-c if(DEBUG)print*,'>>>> try to include a new cluster'
distmin=1000000.
xmm_A = 0. !---------------------------
ymm_A = 0. ! init variables that
@@ -3400,6 +3442,7 @@
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
@@ -3407,14 +3450,20 @@
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,
-c $ 'ETA2','ETA2',
$ PFA,PFA,
- $ AXV_STORE(nplanes-ip+1,ibest),0.)
+ $ AXV_STORE(nplanes-ip+1,ibest),0.,
+ $ bxyz(1),
+ $ bxyz(2)
+ $ )
distance = distance_to(XP,YP)
-c if(DEBUG)print*,'normalized distance ',distance
- if(DEBUG)print*,'( cl-X ',icx
- $ ,' in cp ',id,' ) normalized distance ',distance
+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
@@ -3426,8 +3475,8 @@
rymm = resyPAM
distmin = distance
iclm = icx
-c dedxmm = dedx(icx) !(1)
- dedxmmx = dedx(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
+c dedxmm = sgnl(icx) !(1)
+ dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
dedxmmy = 0. !(1)
endif
11881 continue
@@ -3435,13 +3484,20 @@
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,
-c $ 'ETA2','ETA2',
$ PFA,PFA,
- $ 0.,AYV_STORE(nplanes-ip+1,ibest))
+ $ 0.,AYV_STORE(nplanes-ip+1,ibest),
+ $ bxyz(1),
+ $ bxyz(2)
+ $ )
distance = distance_to(XP,YP)
- if(DEBUG)print*,'( cl-Y ',icy
- $ ,' in cp ',id,' ) normalized distance ',distance
+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
@@ -3453,34 +3509,39 @@
rymm = resyPAM
distmin = distance
iclm = icy
-c dedxmm = dedx(icy) !(1)
+c dedxmm = sgnl(icy) !(1)
dedxmmx = 0. !(1)
- dedxmmy = dedx(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
+ dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
endif
11882 continue
enddo !end loop on cluster inside couples
-*----- single clusters -----------------------------------------------
+*----- single clusters -----------------------------------------------
do ic=1,ncls(ip) !loop on single clusters
icl=cls(ip,ic)
-c if(cl_used(icl).eq.1.or. !if the cluster is already used
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,
-c $ 'ETA2','ETA2',
$ PFA,PFA,
- $ AXV_STORE(nplanes-ip+1,ibest),0.)
+ $ AXV_STORE(nplanes-ip+1,ibest),0.,
+ $ bxyz(1),
+ $ bxyz(2)
+ $ )
else !<---- Y view
call xyz_PAM(0,icl,ist,
-c $ 'ETA2','ETA2',
$ PFA,PFA,
- $ 0.,AYV_STORE(nplanes-ip+1,ibest))
+ $ 0.,AYV_STORE(nplanes-ip+1,ibest),
+ $ bxyz(1),
+ $ bxyz(2)
+ $ )
endif
distance = distance_to(XP,YP)
- if(DEBUG)print*,'( cl-s ',icl
- $ ,' ) normalized distance ',distance
+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
@@ -3492,43 +3553,61 @@
rymm = resyPAM
distmin = distance
iclm = icl
-c dedxmm = dedx(icl) !(1)
if(mod(VIEW(icl),2).eq.0)then !<---- X view
- dedxmmx = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)
- dedxmmy = 0. !(1)
+ dedxmmx = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
+ dedxmmy = 0.
else !<---- Y view
- dedxmmx = 0. !(1)
- dedxmmy = dedx(icl)/mip(VIEW(icl),LADDER(icl)) !(1)(2)
+ dedxmmx = 0.
+ dedxmmy = sgnl(icl)/mip(VIEW(icl),LADDER(icl))
endif
endif
18882 continue
enddo !end loop on single clusters
- if(distmin.le.clinc)then
-
- CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<
-* ----------------------------
+ if(iclm.ne.0)then
if(mod(VIEW(iclm),2).eq.0)then
- XGOOD(nplanes-ip+1)=1.
- resx(nplanes-ip+1)=rxmm
- if(DEBUG)print*,'%%%% included X-cl ',iclm
- $ ,' ( norm.dist.= ',distmin,', cut ',clinc,' )'
- else
- YGOOD(nplanes-ip+1)=1.
- resy(nplanes-ip+1)=rymm
- if(DEBUG)print*,'%%%% included Y-cl ',iclm
- $ ,' ( norm.dist.= ',distmin,', cut ',clinc,' )'
+ 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
-* ----------------------------
- xm_A(nplanes-ip+1) = xmm_A
- ym_A(nplanes-ip+1) = ymm_A
- xm_B(nplanes-ip+1) = xmm_B
- ym_B(nplanes-ip+1) = ymm_B
- zm(nplanes-ip+1) = (zmm_A+zmm_B)/2.
-c dedxtrk(nplanes-ip+1) = dedxmm !<<< !(1)
- dedxtrk_x(nplanes-ip+1) = dedxmmx !<<< !(1)
- dedxtrk_y(nplanes-ip+1) = dedxmmy !<<< !(1)
-* ----------------------------
endif
endif
133 continue
@@ -3539,6 +3618,7 @@
return
end
+
***************************************************
* *
* *
@@ -3547,186 +3627,7 @@
* *
* *
**************************************************
-cccccc 12/08/2006 modified by elena ---> (1)
*
- subroutine clean_XYclouds(ibest,iflag)
-
- include 'commontracker.f'
- include 'level1.f'
- include 'common_momanhough.f'
-c include 'momanhough_init.f'
- include 'level2.f' !(1)
-c include 'calib.f'
-c include 'level1.f'
-
-
-
- do ip=1,nplanes !loop on planes
-
- id=CP_STORE(nplanes-ip+1,ibest)
- icl=CLS_STORE(nplanes-ip+1,ibest)
- if(id.ne.0.or.icl.ne.0)then
- if(id.ne.0)then
- iclx=clx(ip,icp_cp(id))
- icly=cly(ip,icp_cp(id))
-c cl_used(iclx)=1 !tag used clusters
-c cl_used(icly)=1 !tag used clusters
- cl_used(iclx)=ntrk !tag used clusters !(1)
- cl_used(icly)=ntrk !tag used clusters !(1)
- elseif(icl.ne.0)then
-c cl_used(icl)=1 !tag used clusters
- cl_used(icl)=ntrk !tag used clusters !1)
- endif
-
-c if(DEBUG)then
-c print*,ip,' <<< ',id
-c endif
-* -----------------------------
-* remove the couple from clouds
-* remove also vitual couples containing the
-* selected clusters
-* -----------------------------
- do icp=1,ncp_plane(ip)
- if(
- $ clx(ip,icp).eq.iclx
- $ .or.
- $ clx(ip,icp).eq.icl
- $ .or.
- $ cly(ip,icp).eq.icly
- $ .or.
- $ cly(ip,icp).eq.icl
- $ )then
- id=id_cp(ip,icp,1)
- if(DEBUG)then
- print*,ip,' <<< cp ',id
- $ ,' ( cl-x '
- $ ,clx(ip,icp)
- $ ,' cl-y '
- $ ,cly(ip,icp),' ) --> removed'
- endif
-* -----------------------------
-* remove the couple from clouds
- do iyz=1,nclouds_yz
- if(cpcloud_yz(iyz,abs(id)).ne.0)then
- ptcloud_yz(iyz)=ptcloud_yz(iyz)-1
- cpcloud_yz(iyz,abs(id))=0
- endif
- enddo
- do ixz=1,nclouds_xz
- if(cpcloud_xz(ixz,abs(id)).ne.0)then
- ptcloud_xz(ixz)=ptcloud_xz(ixz)-1
- cpcloud_xz(ixz,abs(id))=0
- endif
- enddo
-* -----------------------------
- endif
- enddo
-
- endif
- enddo !end loop on planes
-
- return
- end
-
-
-
-
-c$$$*** * * * *** * * * *** * * * *** * * * *** * * * *** * * * ***
-c$$$ real function fbad_cog(ncog,ic)
-c$$$
-c$$$
-c$$$ include 'commontracker.f'
-c$$$ include 'level1.f'
-c$$$ include 'calib.f'
-c$$$
-c$$$* --> signal of the central strip
-c$$$ sc = CLSIGNAL(INDMAX(ic)) !center
-c$$$
-c$$$* signal of adjacent strips
-c$$$* --> left
-c$$$ sl1 = 0 !left 1
-c$$$ if(
-c$$$ $ (INDMAX(ic)-1).ge.INDSTART(ic)
-c$$$ $ )
-c$$$ $ sl1 = max(0.,CLSIGNAL(INDMAX(ic)-1))
-c$$$
-c$$$ sl2 = 0 !left 2
-c$$$ if(
-c$$$ $ (INDMAX(ic)-2).ge.INDSTART(ic)
-c$$$ $ )
-c$$$ $ sl2 = max(0.,CLSIGNAL(INDMAX(ic)-2))
-c$$$
-c$$$* --> right
-c$$$ sr1 = 0 !right 1
-c$$$ if(
-c$$$ $ (ic.ne.NCLSTR1.and.(INDMAX(ic)+1).lt.INDSTART(ic+1))
-c$$$ $ .or.
-c$$$ $ (ic.eq.NCLSTR1.and.(INDMAX(ic)+1).le.TOTCLLENGTH)
-c$$$ $ )
-c$$$ $ sr1 = max(0.,CLSIGNAL(INDMAX(ic)+1))
-c$$$
-c$$$ sr2 = 0 !right 2
-c$$$ if(
-c$$$ $ (ic.ne.NCLSTR1.and.(INDMAX(ic)+2).lt.INDSTART(ic+1))
-c$$$ $ .or.
-c$$$ $ (ic.eq.NCLSTR1.and.(INDMAX(ic)+2).le.TOTCLLENGTH)
-c$$$ $ )
-c$$$ $ sr2 = max(0.,CLSIGNAL(INDMAX(ic)+2))
-c$$$
-c$$$
-c$$$ if(mod(int(VIEW(ic)),2).eq.1)then !Y-view
-c$$$ f = 4.
-c$$$ si = 8.4
-c$$$ else !X-view
-c$$$ f = 6.
-c$$$ si = 3.9
-c$$$ endif
-c$$$
-c$$$ fbad_cog = 1.
-c$$$ f0 = 1
-c$$$ f1 = 1
-c$$$ f2 = 1
-c$$$ f3 = 1
-c$$$ if(sl1.gt.sr1.and.sl1.gt.0.)then
-c$$$
-c$$$ if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)) ).eq.0)f0=f
-c$$$ if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f1=f
-c$$$c if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f3=f
-c$$$
-c$$$ if(ncog.eq.2.and.sl1.ne.0)then
-c$$$ fbad_cog = (f1**2*sc**2/sl1**2+f0**2)/(sc**2/sl1**2+1.)
-c$$$ elseif(ncog.eq.3.and.sl1.ne.0.and.sr1.ne.0)then
-c$$$ fbad_cog = 1.
-c$$$ elseif(ncog.eq.4.and.sl1.ne.0.and.sr1.ne.0.and.sl2.ne.0)then
-c$$$ fbad_cog = 1.
-c$$$ else
-c$$$ fbad_cog = 1.
-c$$$ endif
-c$$$
-c$$$ elseif(sl1.le.sr1.and.sr1.gt.0.)then
-c$$$
-c$$$
-c$$$ if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)) ).eq.0)f0=f
-c$$$ if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)+1)).eq.0)f1=f
-c$$$c if(BAD(VIEW(ic),nvk(MAXS(ic)),nst(MAXS(ic)-1)).eq.0)f3=f
-c$$$
-c$$$ if(ncog.eq.2.and.sr1.ne.0)then
-c$$$ fbad_cog = (f1**2*sc**2/sr1**2+f0**2)/(sc**2/sr1**2+1.)
-c$$$ elseif(ncog.eq.3.and.sr1.ne.0.and.sl1.ne.0)then
-c$$$ fbad_cog = 1.
-c$$$ elseif(ncog.eq.4.and.sr1.ne.0.and.sl1.ne.0.and.sr2.ne.0)then
-c$$$ fbad_cog = 1.
-c$$$ else
-c$$$ fbad_cog = 1.
-c$$$ endif
-c$$$
-c$$$ endif
-c$$$
-c$$$ fbad_cog = sqrt(fbad_cog)
-c$$$
-c$$$ return
-c$$$ end
-c$$$
@@ -3738,13 +3639,22 @@
include 'level1.f'
include 'common_momanhough.f'
include 'level2.f'
-c include 'level1.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
@@ -3755,12 +3665,23 @@
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
@@ -3780,6 +3701,10 @@
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
@@ -3812,7 +3737,7 @@
alfayz1_nt(idb)=0
alfayz2_nt(idb)=0
enddo
- do itr=1,ntrpl_max_nt
+ do itr=1,ntrpt_max_nt
tr_cloud_nt(itr)=0
alfaxz1_nt(itr)=0
alfaxz2_nt(itr)=0
@@ -3841,7 +3766,7 @@
alfayz1(idb)=0
alfayz2(idb)=0
enddo
- do itr=1,ntrpl_max
+ do itr=1,ntrpt_max
tr_cloud(itr)=0
cpxz1(itr)=0
cpxz2(itr)=0
@@ -3889,19 +3814,27 @@
include 'commontracker.f'
-c include 'level1.f'
include 'level1.f'
include 'common_momanhough.f'
include 'level2.f'
include 'common_mini_2.f'
- real sinth,phi,pig
+ 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 = al(4)
- sinth = al(3)
+* -------------------------------------
+ phi = REAL(al(4))
+ sinth = REAL(al(3))
if(sinth.lt.0)then
sinth = -sinth
phi = phi + pig
@@ -3912,14 +3845,13 @@
$ 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))
@@ -3928,28 +3860,152 @@
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))
- dedx_x(ip,ntr) = sngl(dedxtrk_x(ip)) !(2)
- dedx_y(ip,ntr) = sngl(dedxtrk_y(ip)) !(2)
+ 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)
- id = CP_STORE(ip,IDCAND)
+
+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)
- if(id.ne.0)then
+ 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))
-c print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr)
- elseif(icl.ne.0)then
- if(mod(VIEW(icl),2).eq.0)cltrx(ip,ntr)=icl
- if(mod(VIEW(icl),2).eq.1)cltry(ip,ntr)=icl
-c print*,ip,' ',cltrx(ip,ntr),cltry(ip,ntr)
+
+ 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
@@ -3962,7 +4018,6 @@
* -------------------------------------------------------
include 'commontracker.f'
-c include 'level1.f'
include 'calib.f'
include 'level1.f'
include 'common_momanhough.f'
@@ -3970,54 +4025,70 @@
include 'common_xyzPAM.f'
* count #cluster per plane not associated to any track
-c good2=1!.true.
nclsx = 0
nclsy = 0
do iv = 1,nviews
- if( mask_view(iv).ne.0 )good2(iv) = 20+mask_view(iv)
+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
- ip=nplanes-npl(VIEW(icl))+1
+
if(mod(VIEW(icl),2).eq.0)then !=== X views
+
nclsx = nclsx + 1
planex(nclsx) = ip
- sgnlxs(nclsx) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)
+ 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.)
- call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.)
- xs(is,nclsx) = (xPAM_A+xPAM_B)/2
+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
-c$$$ print*,'nclsx ',nclsx
-c$$$ print*,'planex(nclsx) ',planex(nclsx)
-c$$$ print*,'sgnlxs(nclsx) ',sgnlxs(nclsx)
-c$$$ print*,'xs(1,nclsx) ',xs(1,nclsx)
-c$$$ print*,'xs(2,nclsx) ',xs(2,nclsx)
else !=== Y views
nclsy = nclsy + 1
planey(nclsy) = ip
- sgnlys(nclsy) = dedx(icl)/mip(VIEW(icl),LADDER(icl))!(2)
+ 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.)
- call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.)
- ys(is,nclsy) = (yPAM_A+yPAM_B)/2
+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
-c$$$ print*,'nclsy ',nclsy
-c$$$ print*,'planey(nclsy) ',planey(nclsy)
-c$$$ print*,'sgnlys(nclsy) ',sgnlys(nclsy)
-c$$$ print*,'ys(1,nclsy) ',ys(1,nclsy)
-c$$$ print*,'ys(2,nclsy) ',ys(2,nclsy)
endif
endif
-c print*,icl,cl_used(icl),cl_good(icl),ip,VIEW(icl)!nclsx(ip),nclsy(ip)
***** 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
@@ -4079,7 +4150,7 @@
alfayz2_av_nt(iyz)=alfayz2_av(iyz)
nnn=nnn+ptcloud_yz(iyz)
enddo
- do ipt=1,nnn
+ do ipt=1,min(ndblt_max_nt,nnn)
db_cloud_nt(ipt)=db_cloud(ipt)
enddo
endif
@@ -4092,7 +4163,7 @@
alfaxz3_av_nt(ixz)=alfaxz3_av(ixz)
nnn=nnn+ptcloud_xz(ixz)
enddo
- do ipt=1,nnn
+ do ipt=1,min(ntrpt_max_nt,nnn)
tr_cloud_nt(ipt)=tr_cloud(ipt)
enddo
endif