--- DarthVader/TrackerLevel2/src/F77/analysissubroutines.f 2008/11/25 14:41:38 1.36
+++ DarthVader/TrackerLevel2/src/F77/analysissubroutines.f 2009/12/14 10:51:40 1.43
@@ -20,8 +20,6 @@
include 'calib.f'
include 'level2.f'
-
-
c print*,'======================================================'
c$$$ do ic=1,NCLSTR1
c$$$ if(.false.
@@ -74,12 +72,11 @@
*-------------------------------------------------------------------------------
*-------------------------------------------------------------------------------
-
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
@@ -113,6 +110,7 @@
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
*-------------------------------------------------------------------------------
@@ -163,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(nclouds_yz.eq.0)goto 880 !go to next event
- if(planehit.eq.3) goto 881
+
+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
@@ -243,15 +269,20 @@
*
*-------------------------------------------------------------------------------
*-------------------------------------------------------------------------------
- ntrk=0 !counter of identified physical tracks
+ccc ntrk=0 !counter of identified physical tracks
-11111 continue !<<<<<<< come here when performing a new search
+c11111 continue !<<<<<<< come here when performing a new search
+ continue !<<<<<<< come here when performing a new search
+
+ if(nclouds_xz.eq.0)goto 880 !go to next event
+ if(nclouds_yz.eq.0)goto 880 !go to next event
c iflag=0
call clouds_to_ctrack(iflag)
if(iflag.eq.1)then !no candidate tracks found
goto 880 !fill ntp and go to next event
endif
+ if(ntracks.eq.0)goto 880 !go to next event
FIMAGE=.false. !processing best track (not track image)
ibest=0 !best track among candidates
@@ -272,7 +303,7 @@
* 1st) decreasing n.points
* 2nd) increasing chi**2
* -------------------------------------------------------
- rchi2best=1000000000.
+ rchi2best=1000000000.
ndofbest=0
do i=1,ntracks
ndof=0
@@ -295,24 +326,6 @@
endif
enddo
-c$$$ rchi2best=1000000000.
-c$$$ ndofbest=0 !(1)
-c$$$ do i=1,ntracks
-c$$$ if(RCHI2_STORE(i).lt.rchi2best.and.
-c$$$ $ RCHI2_STORE(i).gt.0)then
-c$$$ ndof=0 !(1)
-c$$$ do ii=1,nplanes !(1)
-c$$$ ndof=ndof !(1)
-c$$$ $ +int(xgood_store(ii,i)) !(1)
-c$$$ $ +int(ygood_store(ii,i)) !(1)
-c$$$ enddo !(1)
-c$$$ if(ndof.ge.ndofbest)then !(1)
-c$$$ ibest=i
-c$$$ rchi2best=RCHI2_STORE(i)
-c$$$ ndofbest=ndof !(1)
-c$$$ endif !(1)
-c$$$ endif
-c$$$ enddo
if(ibest.eq.0)goto 880 !>> no good candidates
*-------------------------------------------------------------------------------
@@ -350,7 +363,6 @@
do i=1,5
AL_GUESS(i)=AL(i)
enddo
-c print*,'## guess: ',al
do i=1,5
AL(i)=dble(AL_STORE(i,icand))
@@ -365,20 +377,14 @@
if(VERBOSE.EQ.1)iprint=1
if(DEBUG.EQ.1)iprint=2
call mini2(jstep,ifail,iprint)
- if(ifail.ne.0) 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 *** (after refinement) '
$ ,iev
-
-c$$$ print*,'guess: ',(al_guess(i),i=1,5)
-c$$$ print*,'previous: ',(al_store(i,icand),i=1,5)
-c$$$ print*,'result: ',(al(i),i=1,5)
-c$$$ print*,'xgood ',xgood
-c$$$ print*,'ygood ',ygood
-c$$$ print*,'----------------------------------------------'
endif
-c chi2=-chi2
endif
if(DEBUG.EQ.1)then
@@ -424,37 +430,14 @@
if(WARNING.EQ.1)print*,'** WARNING ** sensor=0'
goto 122 !jump
endif
-c print*,'is1 ',is1
do ip=1,NPLANES
is2 = SENSOR_STORE(ip,icand) !sensor
-c print*,'is2 ',is2,' ip ',ip
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! '
-* now search for track-image among track candidates
-c$$$ do i=1,ntracks
-c$$$ iimage=i
-c$$$ do ip=1,nplanes
-c$$$ if( CP_STORE(nplanes-ip+1,icand).ne.
-c$$$ $ -1*CP_STORE(nplanes-ip+1,i).and.
-c$$$ $ CP_STORE(nplanes-ip+1,i).ne.0.and.
-c$$$ $ CP_STORE(nplanes-ip+1,icand).ne.0 )iimage=0
-c$$$ print*,' track ',i,' CP ',CP_STORE(nplanes-ip+1,i)
-c$$$ $ ,CP_STORE(nplanes-ip+1,icand),iimage
-c$$$ enddo
-c$$$ if( iimage.ne.0.and.
-c$$$c $ RCHI2_STORE(i).le.CHI2MAX.and.
-c$$$c $ RCHI2_STORE(i).gt.0.and.
-c$$$ $ .true.)then
-c$$$ if(DEBUG.EQ.1)print*,'Track candidate ',iimage
-c$$$ $ ,' >>> TRACK IMAGE >>> of'
-c$$$ $ ,ibest
-c$$$ goto 122 !image track found
-c$$$ endif
-c$$$ enddo
* ---------------------------------------------------------------
* take the candidate with the greatest number of matching couples
* if more than one satisfies the requirement,
@@ -486,8 +469,6 @@
endif
endif
-c$$$ print*,' track ',i,' CP ',CP_STORE(nplanes-ip+1,i)
-c$$$ $ ,CP_STORE(nplanes-ip+1,icand),ncp
enddo
117 continue
trackimage(i)=ncp !number of matching couples
@@ -501,8 +482,9 @@
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.1)then
+ if(ngood.gt.0)then
* -------------------------------------------------------
* order track-candidates according to:
* 1st) decreasing n.points
@@ -532,27 +514,20 @@
endif
endif
enddo
+
+ if(DEBUG.EQ.1)then
+ print*,'Track candidate ',iimage
+ $ ,' >>> TRACK IMAGE >>> of'
+ $ ,ibest
+ endif
endif
- if(DEBUG.EQ.1)print*,'Track candidate ',iimage
- $ ,' >>> TRACK IMAGE >>> of'
- $ ,ibest
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)
-
if(ntrk.eq.NTRKMAX)then
if(verbose.eq.1)
$ print*,
@@ -562,43 +537,29 @@
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.EQ.1)then
- print*,'***** NEW SEARCH ****'
- endif
- goto 11111 !try new search
-
- endif
-* **********************************************
-
-
880 return
end
@@ -699,7 +660,6 @@
parameter (ndivx=30)
-c$$$ print*,icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy
resxPAM = 0
resyPAM = 0
@@ -713,7 +673,6 @@
xPAM_B = 0.D0
yPAM_B = 0.D0
zPAM_B = 0.D0
-cc print*,'## xyz_PAM: ',icx,icy,sensor,PFAx,PFAy,angx,angy
if(sensor.lt.1.or.sensor.gt.2)then
print*,'xyz_PAM ***ERROR*** wrong input '
@@ -757,7 +716,8 @@
stripx = stripx + corr
resxPAM = res
- 10 endif
+ 10 continue
+ endif
* -----------------
* CLUSTER Y
@@ -803,9 +763,9 @@
stripy = stripy + corr
resyPAM = res
- 20 endif
+ 20 continue
+ endif
-cc print*,'## stripx,stripy ',stripx,stripy
c===========================================================
C COUPLE
@@ -885,7 +845,9 @@
nldx = nldy
viewx = viewy + 1
- yi = dcoordsi(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
@@ -895,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===========================================================
@@ -907,8 +867,6 @@
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...
if(DEBUG.EQ.1) then
@@ -916,7 +874,10 @@
$ ' WARNING: false X strip: strip ',stripx
endif
endif
- xi = dcoordsi(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
@@ -932,8 +893,6 @@
yi_B = yi
endif
-c print*,'X-cl ',icx,stripx,' --> ',xi
-c print*,yi_A,' <--> ',yi_B
else
if(DEBUG.EQ.1) then
@@ -979,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
@@ -989,9 +966,12 @@
c in PAMELA reference system
c------------------------------------------------------------------------
- xPAM = 0.D0
- yPAM = 0.D0
- zPAM = 0.D0
+ 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
@@ -1002,7 +982,6 @@
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
if(DEBUG.EQ.1) then
@@ -1013,9 +992,6 @@
endif
-c print*,'## xPAM,yPAM,zPAM ',xPAM,yPAM,zPAM
-c print*,'## xPAM_A,yPAM_A,zPAM_A ',xPAM_A,yPAM_A,zPAM_A
-c print*,'## xPAM_B,yPAM_B,zPAM_B ',xPAM_B,yPAM_B,zPAM_B
100 continue
end
@@ -1063,19 +1039,11 @@
call idtoc(pfaid,PFAx)
call idtoc(pfaid,PFAy)
-cc print*,PFAx,PFAy
-
-c$$$ call xyz_PAM(icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy)
-
-c$$$ print*,icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy
if(icx.ne.0.and.icy.ne.0)then
ipx=npl(VIEW(icx))
ipy=npl(VIEW(icy))
-c$$$ if( (nplanes-ipx+1).ne.ip.or.(nplanes-ipy+1).ne.ip )
-c$$$ $ print*,'xyzpam: ***WARNING*** clusters ',icx,icy
-c$$$ $ ,' does not belong to the correct plane: ',ip,ipx,ipy
if( (nplanes-ipx+1).ne.ip )then
print*,'xyzpam: ***WARNING*** cluster ',icx
@@ -1110,9 +1078,6 @@
elseif(icx.eq.0.and.icy.ne.0)then
ipy=npl(VIEW(icy))
-c$$$ if((nplanes-ipy+1).ne.ip)
-c$$$ $ print*,'xyzpam: ***WARNING*** clusters ',icx,icy
-c$$$ $ ,' does not belong to the correct plane: ',ip,ipx,ipy
if( (nplanes-ipy+1).ne.ip )then
print*,'xyzpam: ***WARNING*** cluster ',icy
$ ,' does not belong to plane: ',ip
@@ -1127,34 +1092,22 @@
resx(ip) = 1000.
resy(ip) = resyPAM
-cPP --- old ---
c$$$ xm(ip) = -100.
c$$$ ym(ip) = -100.
c$$$ zm(ip) = (zPAM_A+zPAM_B)/2.
-c$$$ xm_A(ip) = xPAM_A
-c$$$ ym_A(ip) = yPAM_A
-c$$$ xm_B(ip) = xPAM_B
-c$$$ ym_B(ip) = yPAM_B
-cPP --- new ---
- xm(ip) = -100.
- ym(ip) = -100.
- zm(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
+ xm(ip) = xPAM
+ ym(ip) = yPAM
+ zm(ip) = zPAM
xm_A(ip) = xPAM_A
ym_A(ip) = yPAM_A
- zm_A(ip) = zPAM_A
xm_B(ip) = xPAM_B
ym_B(ip) = yPAM_B
- zm_B(ip) = zPAM_B
-cPP -----------
c zv(ip) = (zPAM_A+zPAM_B)/2.
elseif(icx.ne.0.and.icy.eq.0)then
ipx=npl(VIEW(icx))
-c$$$ if((nplanes-ipx+1).ne.ip)
-c$$$ $ print*,'xyzpam: ***WARNING*** clusters ',icx,icy
-c$$$ $ ,' does not belong to the correct plane: ',ip,ipx,ipy
if( (nplanes-ipx+1).ne.ip )then
print*,'xyzpam: ***WARNING*** cluster ',icx
@@ -1170,26 +1123,17 @@
resx(ip) = resxPAM
resy(ip) = 1000.
-cPP --- old ---
c$$$ xm(ip) = -100.
c$$$ ym(ip) = -100.
c$$$ zm(ip) = (zPAM_A+zPAM_B)/2.
-c$$$ xm_A(ip) = xPAM_A
-c$$$ ym_A(ip) = yPAM_A
-c$$$ xm_B(ip) = xPAM_B
-c$$$ ym_B(ip) = yPAM_B
-cPP --- new ---
- xm(ip) = -100.
- ym(ip) = -100.
- zm(ip) = z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
+ xm(ip) = xPAM
+ ym(ip) = yPAM
+ zm(ip) = zPAM
xm_A(ip) = xPAM_A
ym_A(ip) = yPAM_A
- zm_A(ip) = zPAM_A
xm_B(ip) = xPAM_B
ym_B(ip) = yPAM_B
- zm_B(ip) = zPAM_B
-cPP -----------
-
+
c zv(ip) = (zPAM_A+zPAM_B)/2.
else
@@ -1198,7 +1142,6 @@
if(lad.ne.0)il=lad
is = 1
if(sensor.ne.0)is=sensor
-c print*,nplanes-ip+1,il,is
xgood(ip) = 0.
ygood(ip) = 0.
@@ -1218,12 +1161,10 @@
endif
if(DEBUG.EQ.1)then
-c print*,'----------------------------- track coord'
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)
-c$$$ print*,'-----------------------------'
endif
end
@@ -1266,9 +1207,9 @@
* ----------------------
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.
@@ -1310,16 +1251,13 @@
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.
@@ -1340,17 +1278,9 @@
c$$$ $ ((yPAM-YPP)/resyPAM)**2
distance=dsqrt(distance)
-c$$$ print*,xPAM,yPAM,zPAM
-c$$$ $ ,' --- distance_to --- ',xpp,ypp
-c$$$ print*,' resolution ',resxPAM,resyPAM
else
-c print*
-c $ ,' function distance_to ---> wrong usage!!!'
-c print*,' xPAM,yPAM,zPAM ',xPAM,yPAM,zPAM
-c print*,' xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b '
-c $ ,xPAM_A,yPAM_A,zPAM_A,xPAM_b,yPAM_b,zPAM_b
endif
distance_to = sngl(distance)
@@ -1415,12 +1345,6 @@
c------------------------------------------------------------------------
c (xi,yi,zi) = mechanical coordinates in the silicon sensor frame
c------------------------------------------------------------------------
-c if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
-c $ .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...
-c print*,'whichsensor: ',
-c $ ' WARNING: false X strip: strip ',stripx
-c endif
xi = dcoordsi(stripx,viewx)
yi = dcoordsi(stripy,viewy)
zi = 0.D0
@@ -1448,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.
@@ -1574,7 +1496,6 @@
is_cp=0
if(id.lt.0)is_cp=1
if(id.gt.0)is_cp=2
-c if(id.eq.0)print*,'IS_CP ===> wrong couple id !!!'
return
end
@@ -1673,11 +1594,13 @@
integer iflag
integer badseed,badclx,badcly
-
+
+ iflag = iflag
if(DEBUG.EQ.1)print*,'cl_to_couples:'
+cc if(RECOVER_SINGLETS.and..not.SECOND_SEARCH)goto 80
+
* init variables
- ncp_tot=0
do ip=1,nplanes
do ico=1,ncouplemax
clx(ip,ico)=0
@@ -1690,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
@@ -1719,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)
* ----------------------------------------------------
@@ -1726,7 +1651,7 @@
* ----------------------------------------------------
* cut on charge (X VIEW)
* ----------------------------------------------------
- if(sgnl(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
@@ -1767,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)
@@ -1776,7 +1703,7 @@
* ----------------------------------------------------
* cut on charge (Y VIEW)
* ----------------------------------------------------
- if(sgnl(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
@@ -1854,8 +1781,8 @@
$ ,'( ',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
+ mask_view(nviewx(nplx))= mask_view(nviewx(nplx))+ 2**1
+ mask_view(nviewy(nply))= mask_view(nviewy(nply))+ 2**1
goto 10
endif
@@ -1875,7 +1802,6 @@
10 continue
enddo !end loop on clusters(X)
-
do icl=1,nclstr1
if(cl_single(icl).eq.1)then
ip=npl(VIEW(icl))
@@ -1883,18 +1809,86 @@
cls(ip,ncls(ip))=icl
endif
enddo
+
+c 80 continue
+ continue
if(DEBUG.EQ.1)then
print*,'clusters ',nclstr1
print*,'good ',(cl_good(i),i=1,nclstr1)
+ print*,'used ',(cl_used(i),i=1,nclstr1)
print*,'singlets',(cl_single(i),i=1,nclstr1)
print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)
endif
+
+
+ if(.not.RECOVER_SINGLETS)goto 81
+
+* ////////////////////////////////////////////////
+* PATCH to recover events with less than 3 couples
+* ////////////////////////////////////////////////
+* loop over singlet and create "fake" couples
+* (with clx=0 or cly=0)
+*
+
+ if(DEBUG.EQ.1)
+ $ print*,'>>> Recover singlets '
+ $ ,'(creates fake couples) <<<'
+ do icl=1,nclstr1
+ if(
+ $ cl_single(icl).eq.1.and.
+ $ cl_used(icl).eq.0.and.
+ $ .true.)then
+* ----------------------------------------------------
+* jump masked views
+* ----------------------------------------------------
+ if( mask_view(VIEW(icl)).ne.0 ) goto 21
+ if(mod(VIEW(icl),2).eq.0)then !=== X views
+* ----------------------------------------------------
+* cut on charge (X VIEW)
+* ----------------------------------------------------
+ if(sgnl(icl).lt.dedx_x_min) goto 21
+
+ nplx=npl(VIEW(icl))
+* ------------------> (FAKE) COUPLE <-----------------
+ ncp_plane(nplx) = ncp_plane(nplx) + 1
+ clx(nplx,ncp_plane(nplx))=icl
+ cly(nplx,ncp_plane(nplx))=0
+c$$$ cl_single(icl)=0! I leave the cluster tagged as singlet!!!
+* ----------------------------------------------------
+
+ else !=== Y views
+* ----------------------------------------------------
+* cut on charge (Y VIEW)
+* ----------------------------------------------------
+ if(sgnl(icl).lt.dedx_y_min) goto 21
+
+ nply=npl(VIEW(icl))
+* ------------------> (FAKE) COUPLE <-----------------
+ ncp_plane(nply) = ncp_plane(nply) + 1
+ clx(nply,ncp_plane(nply))=0
+ cly(nply,ncp_plane(nply))=icl
+c$$$ cl_single(icl)=0! I leave the cluster tagged as singlet!!!
+* ----------------------------------------------------
+
+ endif
+ endif !end "single" condition
+ 21 continue
+ enddo !end loop over clusters
+
+ if(DEBUG.EQ.1)
+ $ print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)
+
+
+ 81 continue
- do ip=1,6
+ ncp_tot=0
+ do ip=1,NPLANES
ncp_tot = ncp_tot + ncp_plane(ip)
enddo
+ if(DEBUG.EQ.1)
+ $ print*,'n.couple tot: ',ncp_tot
return
end
@@ -1957,8 +1951,6 @@
* --------------------------------------------
do ip=1,nplanes
if(ncp_plane(ip).gt.ncouplelimit)then
-c mask_view(nviewx(ip)) = 8
-c mask_view(nviewy(ip)) = 8
mask_view(nviewx(ip)) = mask_view(nviewx(ip)) + 2**7
mask_view(nviewy(ip)) = mask_view(nviewy(ip)) + 2**7
endif
@@ -1969,27 +1961,38 @@
ntrpt=0 !number of triplets
do ip1=1,(nplanes-1) !loop on planes - COPPIA 1
+c$$$ print*,'(1) ip ',ip1
+c$$$ $ ,mask_view(nviewx(ip1))
+c$$$ $ ,mask_view(nviewy(ip1))
if( mask_view(nviewx(ip1)).ne.0 .or.
$ mask_view(nviewy(ip1)).ne.0 )goto 10 !skip plane
do is1=1,2 !loop on sensors - COPPIA 1
do icp1=1,ncp_plane(ip1) !loop on COPPIA 1
icx1=clx(ip1,icp1)
icy1=cly(ip1,icp1)
+
+c$$$ print*,'(1) ip ',ip1,' icp ',icp1
+
c call xyz_PAM(icx1,icy1,is1,'COG2','COG2',0.,0.)!(1)
c call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.) !(1)
call xyz_PAM(icx1,icy1,is1,PFAdef,PFAdef,0.,0.,0.,0.)
xm1=xPAM
ym1=yPAM
zm1=zPAM
-c print*,'***',is1,xm1,ym1,zm1
do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2
+c$$$ print*,'(2) ip ',ip2
+c$$$ $ ,mask_view(nviewx(ip2))
+c$$$ $ ,mask_view(nviewy(ip2))
if( mask_view(nviewx(ip2)).ne.0 .or.
$ mask_view(nviewy(ip2)).ne.0 )goto 20 !skip plane
do is2=1,2 !loop on sensors -ndblt COPPIA 2
do icp2=1,ncp_plane(ip2) !loop on COPPIA 2
icx2=clx(ip2,icp2)
icy2=cly(ip2,icp2)
+
+c$$$ print*,'(2) ip ',ip2,' icp ',icp2
+
c call xyz_PAM
c $ (icx2,icy2,is2,'COG2','COG2',0.,0.)!(1)
c call xyz_PAM
@@ -1999,7 +2002,17 @@
xm2=xPAM
ym2=yPAM
zm2=zPAM
-
+
+* ---------------------------------------------------
+* 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)
@@ -2009,15 +2022,19 @@
$ '** 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
-c 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)
@@ -2026,31 +2043,37 @@
alfayz2(ndblt)=(ym1-ym2)/(zm1-zm2)
* y0 (cm)
alfayz1(ndblt)=alfayz2(ndblt)*(zini-zm1)+ym1
-
+
**** -----------------------------------------------****
**** 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(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
@@ -2058,6 +2081,23 @@
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
@@ -2068,83 +2108,131 @@
xm3=xPAM
ym3=yPAM
zm3=zPAM
+
+
* find the circle passing through the three points
-c$$$ call tricircle(3,xp,zp,angp,resp,chi
-c$$$ $ ,xc,zc,radius,iflag)
- iflag = DEBUG
+ 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*ZINI+BX
+
+ endif
+
+ if( .not.SECOND_SEARCH.and.
+ $ radius**2.lt.(ZINI-zc)**2)goto 29
+
* - - - - - - - - - - - - - - - - - - - - - - - - - - - -
* track parameters on X VIEW
* (3 couples needed)
* - - - - - - - - - - - - - - - - - - - - - - - - - - - -
if(ntrpt.eq.ntrpt_max)then
if(verbose.eq.1)print*,
- $ '** warning ** number of identified '//
- $ 'triplets exceeds vector dimention '
- $ ,'( ',ntrpt_max,' )'
-c good2=.false.
-c goto 880 !fill ntp and go to next event
+ $ '** warning **'//
+ $ ' number of identified '//
+ $ 'triplets exceeds'//
+ $ ' vector dimention '
+ $ ,'( ',ntrpt_max,' )'
+c good2=.false.
+c goto 880 !fill ntp and go to next event
do iv=1,nviews
-c mask_view(iv) = 4
- mask_view(iv)=mask_view(iv)+ 2**3
+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-(ZINI-zc)**2)
+ alfaxz2(ntrpt) = (ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)
+ 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-(ZINI-zc)**2)
+ alfaxz2(ntrpt) = -(ZINI-zc)/sqrt(radius**2-(ZINI-zc)**2)
+ 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
- 31 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
@@ -2216,9 +2304,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
@@ -2254,19 +2339,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
@@ -2282,13 +2359,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
@@ -2305,7 +2382,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
@@ -2315,9 +2394,7 @@
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
* ~~~~~~~~~~~~~~~~~
@@ -2349,11 +2426,9 @@
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.EQ.1)then
- print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'
print*,'>>>> cloud ',nclouds_yz,' --- ',npt,' points'
print*,'- alfayz1 ',alfayz1_av(nclouds_yz)
print*,'- alfayz2 ',alfayz2_av(nclouds_yz)
@@ -2362,10 +2437,6 @@
print*,'cpcloud_yz '
$ ,(cpcloud_yz(nclouds_yz,icp),icp=1,ncp_tot)
print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)
-c$$$ print*,'nt-uple: ptcloud_yz(',nclouds_yz,') = '
-c$$$ $ ,ptcloud_yz(nclouds_yz)
-c$$$ print*,'nt-uple: db_cloud(...) = '
-c$$$ $ ,(db_cloud(iii),iii=npt_tot-npt+1,npt_tot)
endif
* >>> NEW CLOUD <<<
* ~~~~~~~~~~~~~~~~~
@@ -2380,9 +2451,7 @@
endif
if(DEBUG.EQ.1)then
- print*,'---------------------- '
print*,'Y-Z total clouds ',nclouds_yz
- print*,' '
endif
@@ -2445,8 +2514,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
@@ -2480,13 +2547,16 @@
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)
-
+
+
* ------------------------------------------------------------------------
* FORCE INCLUSION OF TRIPLETS COMPOSED BY SAME COUPLES, IGNORING THE IMAGE
* ------------------------------------------------------------------------
@@ -2499,7 +2569,6 @@
$ .true.)istrimage=1
if(distance.lt.cutdistxz.or.istrimage.eq.1)then
-c print*,idb1,idb2,distance,' cloud ',nclouds_yz
if(cpxz1(itr2).gt.0)cp_useds2(cpxz1(itr2))=1
if(cpxz1(itr2).lt.0)cp_useds1(-cpxz1(itr2))=1
if(cpxz2(itr2).gt.0)cp_useds2(cpxz2(itr2))=1
@@ -2518,13 +2587,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
@@ -2539,13 +2608,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
@@ -2555,8 +2625,6 @@
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 triplet
* ~~~~~~~~~~~~~~~~~
@@ -2590,8 +2658,7 @@
npt_tot=npt_tot+npt
if(DEBUG.EQ.1)then
- print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'
- print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'
+ print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'
print*,'- alfaxz1 ',alfaxz1_av(nclouds_xz)
print*,'- alfaxz2 ',alfaxz2_av(nclouds_xz)
print*,'- alfaxz3 ',alfaxz3_av(nclouds_xz)
@@ -2600,10 +2667,6 @@
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 <<<
* ~~~~~~~~~~~~~~~~~
@@ -2617,9 +2680,7 @@
endif
if(DEBUG.EQ.1)then
- print*,'---------------------- '
print*,'X-Z total clouds ',nclouds_xz
- print*,' '
endif
@@ -2675,8 +2736,8 @@
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
@@ -2685,28 +2746,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
+ 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
-* get info on
- cpintersec(icp)=min(
- $ cpcloud_yz(iyz,icp),
- $ cpcloud_xz(ixz,icp))
- if(
- $ (cpcloud_yz(iyz,icp).eq.1.and.cpcloud_xz(ixz,icp).eq.2).or.
- $ (cpcloud_yz(iyz,icp).eq.2.and.cpcloud_xz(ixz,icp).eq.1).or.
- $ .false.)cpintersec(icp)=0
+
+ if(.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
@@ -2733,66 +2838,33 @@
do ip=1,nplanes
nplused=nplused+ hit_plane(ip)
enddo
-
-
+
+ if(nplused.lt.3)goto 888 !next combination
+ccc if(nplused.lt.nplxz_min)goto 888 !next combination
+ccc if(nplused.lt.nplyz_min)goto 888 !next combination
+* -----------------------------------------------------------
+* if in RECOVER_SINGLET mode, the two clouds must have
+* at least ONE intersecting real couple
+* -----------------------------------------------------------
+ if(ncp_ok.lt.1)goto 888 !next combination
+
if(DEBUG.EQ.1)then
- print*,'Combination ',iyz,ixz
- $ ,' db ',ptcloud_yz(iyz)
- $ ,' tr ',ptcloud_xz(ixz)
- $ ,' -----> # matching couples ',ncp_ok
- endif
-
-c if(nplused.lt.nplxz_min)goto 888 !next combination
- if(nplused.lt.nplyz_min)goto 888 !next combination
- if(ncp_ok.lt.ncpok)goto 888 !next combination
-
-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 <-------
-cccc SBAGLIATO
-c$$$ AL_INI(1) = dreal(alfaxz1_av(ixz))
-c$$$ AL_INI(2) = dreal(alfayz1_av(iyz))
-c$$$ AL_INI(4) = PIGR + datan(dreal(alfayz2_av(iyz))
-c$$$ $ /dreal(alfaxz2_av(ixz)))
-c$$$ tath = -dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))
-c$$$ AL_INI(3) = tath/sqrt(1+tath**2)
-c$$$ AL_INI(5) = (1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.
-cccc GIUSTO (ma si sua guess())
-c$$$ AL_INI(1) = dreal(alfaxz1_av(ixz))
-c$$$ AL_INI(2) = dreal(alfayz1_av(iyz))
-c$$$ tath = -dreal(alfaxz2_av(ixz))/dcos(AL_INI(4))
-c$$$ AL_INI(3) = tath/sqrt(1+tath**2)
-c$$$ IF(alfaxz2_av(ixz).NE.0)THEN
-c$$$ AL_INI(4) = PIGR + datan(dreal(alfayz2_av(iyz))
-c$$$ $ /dreal(alfaxz2_av(ixz)))
-c$$$ ELSE
-c$$$ AL_INI(4) = acos(-1.)/2
-c$$$ IF(alfayz2_av(iyz).LT.0)AL_INI(4) = AL_INI(4)+acos(-1.)
-c$$$ ENDIF
-c$$$ IF(alfaxz2_av(ixz).LT.0)AL_INI(4)= acos(-1.)+ AL_INI(4)
-c$$$ AL_INI(4) = -acos(-1.) + AL_INI(4) !from incidence direction to tracking rs
-c$$$
-c$$$ AL_INI(5) = (1.e2*alfaxz3_av(ixz))/(0.3*0.43) !0.
-c$$$
-c$$$ if(AL_INI(5).gt.defmax)goto 888 !next cloud
+ 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*,'track candidate', ntracks+1
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))
@@ -2825,6 +2897,16 @@
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
@@ -2873,13 +2955,42 @@
$ 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
@@ -2920,6 +3031,8 @@
* **********************************************************
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
@@ -2946,7 +3059,7 @@
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))
@@ -2965,11 +3078,19 @@
$ cp_match(ip,hit_plane(ip))
SENSOR_STORE(nplanes-ip+1,ntracks)
$ = is_cp(cp_match(ip,hit_plane(ip)))
- LADDER_STORE(nplanes-ip+1,ntracks)
- $ = LADDER(
+
+ 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
@@ -3003,18 +3124,9 @@
if(ntracks.eq.0)then
iflag=1
- return
+cc return
endif
-c$$$ if(DEBUG.EQ.1)then
-c$$$ print*,'****** TRACK CANDIDATES ***********'
-c$$$ print*,'# R. chi2 RIG'
-c$$$ do i=1,ntracks
-c$$$ print*,i,' --- ',rchi2_store(i),' --- '
-c$$$ $ ,1./abs(AL_STORE(5,i))
-c$$$ enddo
-c$$$ print*,'***********************************'
-c$$$ endif
if(DEBUG.EQ.1)then
print*,'****** TRACK CANDIDATES *****************'
print*,'# R. chi2 RIG ndof'
@@ -3094,7 +3206,42 @@
* 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)
@@ -3119,16 +3266,48 @@
$ bxyz(2)
$ )
- 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
-
- dedxtrk_x(nplanes-ip+1)=sgnl(icx)/mip(VIEW(icx),LADDER(icx))
- dedxtrk_y(nplanes-ip+1)=sgnl(icy)/mip(VIEW(icy),LADDER(icy))
+ 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
* |||||||||||||||||||||||||||||||||||||||||||||||||
* -------------------------------------------------
@@ -3140,6 +3319,10 @@
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
@@ -3162,8 +3345,7 @@
* ===========================================
* STEP 1 >>>>>>> try to include a new couple
* ===========================================
-c if(DEBUG.EQ.1)print*,'>>>> try to include a new couple'
- distmin=1000000.
+ distmin=100000000.
xmm = 0.
ymm = 0.
zmm = 0.
@@ -3176,6 +3358,7 @@
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
@@ -3194,7 +3377,8 @@
distance = distance_to(XP,YP)
c distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
id=id_cp(ip,icp,ist)
- if(DEBUG.EQ.1)print*,'( couple ',id
+ if(DEBUG.EQ.1)
+ $ print*,'( couple ',id
$ ,' ) distance ',distance
if(distance.lt.distmin)then
xmm = xPAM
@@ -3206,14 +3390,12 @@
idm = id
dedxmmx = sgnl(icx)/mip(VIEW(icx),LADDER(icx)) !(1)(2)
dedxmmy = sgnl(icy)/mip(VIEW(icy),LADDER(icy)) !(1)(2)
-c QUIQUI --> non devo moltiplicare per clinc?!?!?!
- clincnewc=10*sqrt(rymm**2+rxmm**2 !QUIQUI
- $ +RCHI2_STORE(ibest)*k(ip)*(cov(1,1)+cov(2,2))) !QUIQUI
+ clincnewc=10*sqrt(rymm**2+rxmm**2
+ $ +RCHI2_STORE(ibest)*k(ip)*(cov(1,1)+cov(2,2)))
endif
1188 continue
enddo !end loop on couples on plane icp
-c if(distmin.le.clinc)then !QUIQUI
- if(distmin.le.clincnewc)then !QUIQUI
+ if(distmin.le.clincnewc)then
* -----------------------------------
xm(nplanes-ip+1) = xmm !<<<
ym(nplanes-ip+1) = ymm !<<<
@@ -3227,14 +3409,13 @@
* -----------------------------------
CP_STORE(nplanes-ip+1,ibest)=idm
if(DEBUG.EQ.1)print*,'%%%% included couple ',idm
- $ ,' (dist.= ',distmin,', cut ',clinc,' )'
+ $ ,' (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.EQ.1)print*,'>>>> try to include a new cluster'
distmin=1000000.
xmm_A = 0. !---------------------------
ymm_A = 0. ! init variables that
@@ -3253,6 +3434,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
@@ -3271,7 +3453,8 @@
$ )
distance = distance_to(XP,YP)
c distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
- if(DEBUG.EQ.1)print*,'( cl-X ',icx
+ if(DEBUG.EQ.1)
+ $ print*,'( cl-X ',icx
$ ,' in cp ',id,' ) distance ',distance
if(distance.lt.distmin)then
xmm_A = xPAM_A
@@ -3304,7 +3487,8 @@
$ )
distance = distance_to(XP,YP)
c distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
- if(DEBUG.EQ.1)print*,'( cl-Y ',icy
+ if(DEBUG.EQ.1)
+ $ print*,'( cl-Y ',icy
$ ,' in cp ',id,' ) distance ',distance
if(distance.lt.distmin)then
xmm_A = xPAM_A
@@ -3324,11 +3508,8 @@
11882 continue
enddo !end loop on cluster inside couples
*----- single clusters -----------------------------------------------
-c print*,'## ncls(',ip,') ',ncls(ip)
do ic=1,ncls(ip) !loop on single clusters
-c print*,'-',ic,'-'
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
@@ -3350,10 +3531,10 @@
distance = distance_to(XP,YP)
c distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
- if(DEBUG.EQ.1)print*,'( cl-s ',icl
+ if(DEBUG.EQ.1)
+ $ print*,'( cl-s ',icl
$ ,' ) distance ',distance
if(distance.lt.distmin)then
-c if(DEBUG.EQ.1)print*,'YES'
xmm_A = xPAM_A
ymm_A = yPAM_A
zmm_A = zPAM_A
@@ -3374,10 +3555,7 @@
endif
18882 continue
enddo !end loop on single clusters
-c print*,'## distmin ', distmin,' clinc ',clinc
-c QUIQUI------------
-c anche qui: non ci vuole clinc???
if(iclm.ne.0)then
if(mod(VIEW(iclm),2).eq.0)then
clincnew=
@@ -3388,40 +3566,34 @@
$ 10*
$ sqrt(rymm**2+RCHI2_STORE(ibest)*k(ip)*cov(2,2))
endif
-c QUIQUI------------
-
- if(distmin.le.clincnew)then !QUIQUI
-c if(distmin.le.clinc)then !QUIQUI
+
+ if(distmin.le.clincnew)then
CLS_STORE(nplanes-ip+1,ibest)=iclm !<<<<
* ----------------------------
-c print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
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
+ if(DEBUG.EQ.1)
+ $ print*,'%%%% included X-cl ',iclm
$ ,'( chi^2, ',RCHI2_STORE(ibest)
$ ,', dist.= ',distmin
- $ ,', cut ',clinc,' )'
+ $ ,', cut ',clincnew,' )'
else
YGOOD(nplanes-ip+1)=1.
resy(nplanes-ip+1)=rymm
- if(DEBUG.EQ.1)print*,'%%%% included Y-cl ',iclm
+ if(DEBUG.EQ.1)
+ $ print*,'%%%% included Y-cl ',iclm
$ ,'( chi^2, ',RCHI2_STORE(ibest)
$ ,', dist.= ', distmin
- $ ,', cut ',clinc,' )'
+ $ ,', cut ',clincnew,' )'
endif
-c print*,'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
* ----------------------------
xm_A(nplanes-ip+1) = xmm_A
ym_A(nplanes-ip+1) = ymm_A
- zm_A(nplanes-ip+1) = zmm_A
xm_B(nplanes-ip+1) = xmm_B
ym_B(nplanes-ip+1) = ymm_B
- zm_B(nplanes-ip+1) = zmm_B
zm(nplanes-ip+1) = (zmm_A+zmm_B)/2.
-c$$$ zm(nplanes-ip+1) =
-c$$$ $ z_mech_sensor(nplanes-ip+1,il,is)*1000./1.d4
dedxtrk_x(nplanes-ip+1) = dedxmmx !<<<
dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
* ----------------------------
@@ -3436,6 +3608,7 @@
return
end
+
***************************************************
* *
* *
@@ -3445,78 +3618,6 @@
* *
**************************************************
*
- subroutine clean_XYclouds(ibest,iflag)
-
- include 'commontracker.f'
- include 'level1.f'
- include 'common_momanhough.f'
- include 'level2.f'
-
- if(DEBUG.EQ.1)print*,'clean_XYclouds:'
-
- 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)=ntrk !tag used clusters
-c$$$ cl_used(icly)=ntrk !tag used clusters
- elseif(icl.ne.0)then
-c$$$ cl_used(icl)=ntrk !tag used clusters
- 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.EQ.1)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
-
-
-
@@ -3716,6 +3817,8 @@
integer ssensor,sladder
pig=acos(-1.)
+
+
* -------------------------------------
chi2_nt(ntr) = sngl(chi2)
nstep_nt(ntr) = nstep
@@ -3763,6 +3866,9 @@
dedx_x(ip,ntr) = sngl(dedxtrk_x(ip)/factor)
dedx_y(ip,ntr) = sngl(dedxtrk_y(ip)/factor)
+
+ccc print*,ip,'dedx >>> ',dedx_x(ip,ntr),dedx_y(ip,ntr)
+
ax = axv_nt(ip,ntr)
ay = ayv_nt(ip,ntr)
bfx = BX_STORE(ip,IDCAND)
@@ -3778,7 +3884,6 @@
angy = effectiveangle(ay,2*ip-1,bfx)
-c print*,'* ',ip,bfx,bfy,angx,angy
id = CP_STORE(ip,IDCAND) ! couple id
icl = CLS_STORE(ip,IDCAND)
@@ -3789,39 +3894,59 @@
if(ip.eq.6.and.ssensor.ne.0)ssensor = 3 - ssensor !notazione paolo x align
LS(IP,ntr) = ssensor+10*sladder
- if(id.ne.0)then
+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
+ if(
+ $ xgood_nt(ip,ntr).eq.1.and.ygood_nt(ip,ntr).eq.1
+ $ .and.
+ $ id.ne.0)then
+
c >>> is a couple
cltrx(ip,ntr) = clx(nplanes-ip+1,icp_cp(id))
cltry(ip,ntr) = cly(nplanes-ip+1,icp_cp(id))
- cl_used(cltrx(ip,ntr)) = 1 !tag used clusters
- cl_used(cltry(ip,ntr)) = 1 !tag used clusters
+ if(clx(nplanes-ip+1,icp_cp(id)).gt.0)then
- xbad(ip,ntr)= nbadstrips(4,clx(nplanes-ip+1,icp_cp(id)))
- ybad(ip,ntr)= nbadstrips(4,cly(nplanes-ip+1,icp_cp(id)))
+ 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)
- if(nsatstrips(cly(nplanes-ip+1,icp_cp(id))).gt.0)
- $ dedx_y(ip,ntr)=-dedx_y(ip,ntr)
+ 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
- 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
- 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
+ cl_used(cltry(ip,ntr)) = 1 !tag used clusters
- elseif(icl.ne.0)then
+ 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
+
+ elseif(icl.ne.0)then
cl_used(icl) = 1 !tag used clusters
@@ -3863,15 +3988,12 @@
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
-c$$$ print*,(xgood(i),i=1,6)
-c$$$ print*,(ygood(i),i=1,6)
-c$$$ print*,(ls(i,ntr),i=1,6)
-c$$$ print*,(dedx_x(i,ntr),i=1,6)
-c$$$ print*,(dedx_y(i,ntr),i=1,6)
-c$$$ print*,'-----------------------'
-
end
subroutine fill_level2_siglets
@@ -3908,8 +4030,6 @@
if(cl_used(icl).eq.0)then !cluster not included in any track
-cc print*,' ic ',icl,' --- stored '
-
if(mod(VIEW(icl),2).eq.0)then !=== X views
nclsx = nclsx + 1
@@ -3920,7 +4040,6 @@
sxbad(nclsx) = nbadstrips(1,icl)
multmaxsx(nclsx) = maxs(icl)+10000*mult(icl)
-cc print*,icl,' >>>> ',sxbad(nclsx)
do is=1,2
c call xyz_PAM(icl,0,is,'COG1',' ',0.,0.)
@@ -3928,11 +4047,6 @@
call xyz_PAM(icl,0,is,PFAdef,' ',0.,0.,0.,0.)
xs(is,nclsx) = (xPAM_A+xPAM_B)/2
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
@@ -3942,7 +4056,6 @@
sybad(nclsy) = nbadstrips(1,icl)
multmaxsy(nclsy) = maxs(icl)+10000*mult(icl)
-cc print*,icl,' >>>> ',sybad(nclsy)
do is=1,2
c call xyz_PAM(0,icl,is,' ','COG1',0.,0.)
@@ -3950,11 +4063,6 @@
call xyz_PAM(0,icl,is,' ',PFAdef,0.,0.,0.,0.)
ys(is,nclsy) = (yPAM_A+yPAM_B)/2
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
@@ -4029,7 +4137,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
@@ -4042,7 +4150,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