--- DarthVader/TrackerLevel2/src/F77/analysissubroutines.f 2007/05/09 07:50:58 1.22
+++ DarthVader/TrackerLevel2/src/F77/analysissubroutines.f 2007/08/28 13:25:50 1.30
@@ -215,6 +215,7 @@
c include 'momanhough_init.f'
logical FIMAGE !
+ real trackimage(NTRACKSMAX)
real*8 AL_GUESS(5)
*-------------------------------------------------------------------------------
@@ -331,7 +332,7 @@
iimage=0
endif
if(icand.eq.0)then
- if(VERBOSE)then
+ if(VERBOSE.EQ.1)then
print*,'HAI FATTO UN CASINO!!!!!! icand = ',icand
$ ,ibest,iimage
endif
@@ -360,12 +361,12 @@
jstep=0 !# minimization steps
iprint=0
-c if(DEBUG)iprint=1
- if(VERBOSE)iprint=1
- if(DEBUG)iprint=2
+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(VERBOSE)then
+ if(VERBOSE.EQ.1)then
print *,
$ '*** MINIMIZATION FAILURE *** (after refinement) '
$ ,iev
@@ -380,7 +381,7 @@
c 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
@@ -391,7 +392,7 @@
endif
c rchi2=chi2/dble(ndof)
- if(DEBUG)then
+ if(DEBUG.EQ.1)then
print*,' '
print*,'****** SELECTED TRACK *************'
print*,'# R. chi2 RIG'
@@ -407,25 +408,140 @@
* ------------- 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
+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,
+* 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
+c$$$ print*,' track ',i,' CP ',CP_STORE(nplanes-ip+1,i)
+c$$$ $ ,CP_STORE(nplanes-ip+1,icand),ncp
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
- $ ,' >>> TRACK IMAGE >>> of'
- $ ,ibest
- goto 122 !image track found
+ 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 there are, choose the best one
+ if(ngood.gt.1)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
+
+ 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
@@ -438,7 +554,7 @@
c $ ,iimage,fimage,ntrk,image(ntrk)
if(ntrk.eq.NTRKMAX)then
- if(verbose)
+ if(verbose.eq.1)
$ print*,
$ '** warning ** number of identified '//
$ 'tracks exceeds vector dimension '
@@ -474,7 +590,7 @@
$ rchi2best.le.CHI2MAX.and.
c $ rchi2best.lt.15..and.
$ .true.)then
- if(DEBUG)then
+ if(DEBUG.EQ.1)then
print*,'***** NEW SEARCH ****'
endif
goto 11111 !try new search
@@ -578,6 +694,9 @@
parameter (ndivx=30)
+
+
+c$$$ print*,icx,icy,sensor,PFAx,PFAy,ax,ay,bfx,bfy
resxPAM = 0
resyPAM = 0
@@ -593,10 +712,16 @@
zPAM_B = 0.
c 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 '
+ print*,'sensor ',sensor
+ icx=0
+ icy=0
+ endif
+
* -----------------
* CLUSTER X
-* -----------------
-
+* -----------------
if(icx.ne.0)then
viewx = VIEW(icx)
@@ -604,6 +729,21 @@
nplx = npl(VIEW(icx))
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
* --------------------------
@@ -664,34 +804,43 @@
elseif(PFAx.eq.'ETA2')then
stripx = stripx + pfaeta2(icx,angx)
- resxPAM = risx_eta2(abs(angx))
+ resxPAM = risxeta2(abs(angx))
resxPAM = resxPAM*fbad_cog(2,icx)
- if(DEBUG.and.fbad_cog(2,icx).ne.1)
- $ print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
+c$$$ if(DEBUG.EQ.1.and.fbad_cog(2,icx).ne.1)
+c$$$ $ print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
elseif(PFAx.eq.'ETA3')then
stripx = stripx + pfaeta3(icx,angx)
- resxPAM = risx_eta3(abs(angx))
+ resxPAM = risxeta3(abs(angx))
resxPAM = resxPAM*fbad_cog(3,icx)
- if(DEBUG.and.fbad_cog(3,icx).ne.1)
- $ print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)
+c if(DEBUG.and.fbad_cog(3,icx).ne.1)
+c $ print*,'BAD icx >>> ',viewx,fbad_cog(3,icx)
elseif(PFAx.eq.'ETA4')then
stripx = stripx + pfaeta4(icx,angx)
- resxPAM = risx_eta4(abs(angx))
+ resxPAM = risxeta4(abs(angx))
resxPAM = resxPAM*fbad_cog(4,icx)
- if(DEBUG.and.fbad_cog(4,icx).ne.1)
- $ print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)
+c if(DEBUG.and.fbad_cog(4,icx).ne.1)
+c $ print*,'BAD icx >>> ',viewx,fbad_cog(4,icx)
elseif(PFAx.eq.'ETA')then
stripx = stripx + pfaeta(icx,angx)
- resxPAM = ris_eta(icx,angx)
+c resxPAM = riseta(icx,angx)
+ resxPAM = riseta(viewx,angx)
resxPAM = resxPAM*fbad_eta(icx,angx)
- if(DEBUG.and.fbad_cog(2,icx).ne.1)
- $ print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
+c if(DEBUG.and.fbad_cog(2,icx).ne.1)
+c $ print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
+
+ elseif(PFAx.eq.'ETAL')then
+
+ stripx = stripx + pfaetal(icx,angx)
+ resxPAM = riseta(viewx,angx)
+ resxPAM = resxPAM*fbad_eta(icx,angx)
+c if(DEBUG.and.fbad_cog(2,icx).ne.1)
+c $ print*,'BAD icx >>> ',viewx,fbad_cog(2,icx)
elseif(PFAx.eq.'COG')then
@@ -700,7 +849,7 @@
resxPAM = resxPAM*fbad_cog(0,icx)
else
- if(DEBUG) print*,'*** Non valid p.f.a. (x) --> ',PFAx
+ if(DEBUG.EQ.1) print*,'*** Non valid p.f.a. (x) --> ',PFAx
endif
@@ -710,13 +859,14 @@
if( nsatstrips(icx).gt.0 )then
stripx = stripx + cog(4,icx)
resxPAM = pitchX*1e-4/sqrt(12.)
- cc=cog(4,icx)
+cc cc=cog(4,icx)
c$$$ print*,icx,' *** ',cc
c$$$ print*,icx,' *** ',resxPAM
endif
- endif
-
+ 10 endif
+
+
* -----------------
* CLUSTER Y
* -----------------
@@ -729,8 +879,22 @@
resyPAM = RESYAV
stripy = float(MAXS(icy))
+ if(
+ $ viewy.lt.1.or.
+ $ viewy.gt.12.or.
+ $ nldy.lt.1.or.
+ $ nldy.gt.3.or.
+ $ stripy.lt.1.or.
+ $ stripy.gt.3072.or.
+ $ .false.)then
+ print*,'xyz_PAM ***ERROR*** wrong input '
+ print*,'icy ',icy,'view ',viewy,'nld ',nldy,'strip ',stripy
+ icy = 0
+ goto 20
+ endif
+
if(icx.ne.0.and.(nply.ne.nplx.or.nldy.ne.nldx))then
- if(DEBUG) then
+ if(DEBUG.EQ.1) then
print*,'xyz_PAM ***ERROR*** invalid cluster couple!!! '
$ ,icx,icy
endif
@@ -782,32 +946,41 @@
elseif(PFAy.eq.'ETA2')then
stripy = stripy + pfaeta2(icy,angy)
- resyPAM = risy_eta2(abs(angy))
+ resyPAM = risyeta2(abs(angy))
resyPAM = resyPAM*fbad_cog(2,icy)
- if(DEBUG.and.fbad_cog(2,icy).ne.1)
- $ print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
+c if(DEBUG.and.fbad_cog(2,icy).ne.1)
+c $ print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
elseif(PFAy.eq.'ETA3')then
stripy = stripy + pfaeta3(icy,angy)
resyPAM = resyPAM*fbad_cog(3,icy)
- if(DEBUG.and.fbad_cog(3,icy).ne.1)
- $ print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)
+c if(DEBUG.and.fbad_cog(3,icy).ne.1)
+c $ print*,'BAD icy >>> ',viewy,fbad_cog(3,icy)
elseif(PFAy.eq.'ETA4')then
stripy = stripy + pfaeta4(icy,angy)
resyPAM = resyPAM*fbad_cog(4,icy)
- if(DEBUG.and.fbad_cog(4,icy).ne.1)
- $ print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)
+c if(DEBUG.and.fbad_cog(4,icy).ne.1)
+c $ print*,'BAD icy >>> ',viewy,fbad_cog(4,icy)
elseif(PFAy.eq.'ETA')then
stripy = stripy + pfaeta(icy,angy)
- resyPAM = ris_eta(icy,angy)
+c resyPAM = riseta(icy,angy)
+ resyPAM = riseta(viewy,angy)
resyPAM = resyPAM*fbad_eta(icy,angy)
- if(DEBUG.and.fbad_cog(2,icy).ne.1)
- $ print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
+c if(DEBUG.and.fbad_cog(2,icy).ne.1)
+c $ print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
+
+ elseif(PFAy.eq.'ETAL')then
+
+ stripy = stripy + pfaetal(icy,angy)
+ resyPAM = riseta(viewy,angy)
+ resyPAM = resyPAM*fbad_eta(icy,angy)
+c if(DEBUG.and.fbad_cog(2,icy).ne.1)
+c $ print*,'BAD icy >>> ',viewy,fbad_cog(2,icy)
elseif(PFAy.eq.'COG')then
@@ -816,7 +989,7 @@
resyPAM = resyPAM*fbad_cog(0,icy)
else
- if(DEBUG) print*,'*** Non valid p.f.a. (x) --> ',PFAx
+ if(DEBUG.EQ.1) print*,'*** Non valid p.f.a. (x) --> ',PFAx
endif
@@ -826,15 +999,15 @@
if( nsatstrips(icy).gt.0 )then
stripy = stripy + cog(4,icy)
resyPAM = pitchY*1e-4/sqrt(12.)
- cc=cog(4,icy)
+cc cc=cog(4,icy)
c$$$ print*,icy,' *** ',cc
c$$$ print*,icy,' *** ',resyPAM
endif
- endif
+ 20 endif
-c print*,'## stripx,stripy ',stripx,stripy
+c$$$ print*,'## stripx,stripy ',stripx,stripy
c===========================================================
C COUPLE
@@ -846,7 +1019,7 @@
c------------------------------------------------------------------------
if(((mod(int(stripx+0.5)-1,1024)+1).le.3)
$ .or.((mod(int(stripx+0.5)-1,1024)+1).ge.1022)) then !X has 1018 strips...
- if(DEBUG) then
+ if(DEBUG.EQ.1) then
print*,'xyz_PAM (couple):',
$ ' WARNING: false X strip: strip ',stripx
endif
@@ -941,7 +1114,7 @@
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) then
+ if(DEBUG.EQ.1) then
print*,'xyz_PAM (X-singlet):',
$ ' WARNING: false X strip: strip ',stripx
endif
@@ -966,7 +1139,7 @@
c print*,yi_A,' <--> ',yi_B
else
- if(DEBUG) then
+ if(DEBUG.EQ.1) then
print *,'routine xyz_PAM ---> not properly used !!!'
print *,'icx = ',icx
print *,'icy = ',icy
@@ -1035,7 +1208,7 @@
c print*,'A-(',xPAM_A,yPAM_A,') B-(',xPAM_B,yPAM_B,')'
else
- if(DEBUG) then
+ if(DEBUG.EQ.1) then
print *,'routine xyz_PAM ---> not properly used !!!'
print *,'icx = ',icx
print *,'icy = ',icy
@@ -1050,6 +1223,188 @@
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
+ $ ,' does not exists (nclstr1=',nclstr1,')'
+ icx = -1*icx
+ icy = -1*icy
+ return
+
+ endif
+
+ 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
+ $ ,' 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.
+ ym_A(ip) = 0.
+ xm_B(ip) = 0.
+ ym_B(ip) = 0.
+
+c zv(ip) = zPAM
+
+ 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
+ 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
+
+ xm(ip) = -100.
+ ym(ip) = -100.
+ zm(ip) = (zPAM_A+zPAM_B)/2.
+ 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))
+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
+ $ ,' 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.
+
+ xm(ip) = -100.
+ ym(ip) = -100.
+ zm(ip) = (zPAM_A+zPAM_B)/2.
+ 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
+c print*,nplanes-ip+1,il,is
+
+ 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
+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
********************************************************************************
********************************************************************************
@@ -1492,6 +1847,8 @@
integer iflag
integer badseed,badclx,badcly
+
+ if(DEBUG.EQ.1)print*,'cl_to_couples:'
* init variables
ncp_tot=0
@@ -1522,9 +1879,11 @@
* mask views with too many clusters
do iv=1,nviews
if( ncl_view(iv).gt. nclusterlimit)then
- mask_view(iv) = 1
- if(DEBUG)print*,' * WARNING * cl_to_couple: n.clusters > '
- $ ,nclusterlimit,' on view ', iv,' --> masked!'
+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
@@ -1662,13 +2021,15 @@
endif
if(ncp_plane(nplx).gt.ncouplemax)then
- if(verbose)print*,
+ if(verbose.eq.1)print*,
$ '** warning ** number of identified '//
$ 'couples on plane ',nplx,
$ 'exceeds vector dimention '
$ ,'( ',ncouplemax,' ) --> masked!'
- mask_view(nviewx(nplx)) = 2
- mask_view(nviewy(nply)) = 2
+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
@@ -1698,10 +2059,10 @@
enddo
- 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*,'singlets',(cl_single(i),i=1,nclstr1)
print*,'couples per plane: ',(ncp_plane(ip),ip=1,nplanes)
endif
@@ -1761,6 +2122,7 @@
real xc,zc,radius
* -----------------------------
+ if(DEBUG.EQ.1)print*,'cp_to_doubtrip:'
* --------------------------------------------
* put a limit to the maximum number of couples
@@ -1769,8 +2131,10 @@
* --------------------------------------------
do ip=1,nplanes
if(ncp_plane(ip).gt.ncouplelimit)then
- mask_view(nviewx(ip)) = 8
- mask_view(nviewy(ip)) = 8
+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
enddo
@@ -1796,8 +2160,7 @@
do ip2=(ip1+1),nplanes !loop on planes - COPPIA 2
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 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)
@@ -1816,14 +2179,15 @@
* (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
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
@@ -1879,6 +2243,9 @@
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
call tricircle(3,xp,zp,angp,resp,chi
$ ,xc,zc,radius,iflag)
c print*,xc,zc,radius
@@ -1895,14 +2262,15 @@
* (3 couples needed)
* - - - - - - - - - - - - - - - - - - - - - - - - - - - -
if(ntrpt.eq.ntrpt_max)then
- if(verbose)print*,
+ 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
@@ -1956,7 +2324,7 @@
10 continue
enddo !end loop on planes - COPPIA 1
- if(DEBUG)then
+ if(DEBUG.EQ.1)then
print*,'--- doublets ',ndblt
print*,'--- triplets ',ntrpt
endif
@@ -2003,6 +2371,7 @@
integer cp_useds1(ncouplemaxtot) ! sensor 1
integer cp_useds2(ncouplemaxtot) ! sensor 2
+ if(DEBUG.EQ.1)print*,'doub_to_YZcloud:'
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
* classification of DOUBLETS
@@ -2129,14 +2498,15 @@
* >>> 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
@@ -2156,14 +2526,16 @@
c print*,'>> ',ipt,db_all(ipt)
enddo
npt_tot=npt_tot+npt
- if(DEBUG)then
+ if(DEBUG.EQ.1)then
print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'
print*,'>>>> cloud ',nclouds_yz,' --- ',npt,' points'
- print*,'- alfayz1 ',alfayz1_av(nclouds_yz)
- print*,'- alfayz2 ',alfayz2_av(nclouds_yz)
- print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot)
- print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot)
- print*,'hit_plane ',(hit_plane(ip),ip=1,nplanes)
+ 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)
c$$$ print*,'nt-uple: ptcloud_yz(',nclouds_yz,') = '
c$$$ $ ,ptcloud_yz(nclouds_yz)
c$$$ print*,'nt-uple: db_cloud(...) = '
@@ -2181,7 +2553,7 @@
goto 90
endif
- if(DEBUG)then
+ if(DEBUG.EQ.1)then
print*,'---------------------- '
print*,'Y-Z total clouds ',nclouds_yz
print*,' '
@@ -2230,6 +2602,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
@@ -2287,7 +2661,18 @@
$ +((alfaxz2(itrref)-alfaxz2(itr2))/Dalfaxz2)**2
distance = sqrt(distance)
- if(distance.lt.cutdistxz)then
+* ------------------------------------------------------------------------
+* 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
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
@@ -2351,14 +2736,15 @@
* ~~~~~~~~~~~~~~~~~
* >>> 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
@@ -2377,14 +2763,16 @@
enddo
npt_tot=npt_tot+npt
- if(DEBUG)then
+ if(DEBUG.EQ.1)then
print*,'-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~'
print*,'>>>> cloud ',nclouds_xz,' --- ',npt,' points'
- print*,'- alfaxz1 ',alfaxz1_av(nclouds_xz)
- print*,'- alfaxz2 ',alfaxz2_av(nclouds_xz)
- print*,'- alfaxz3 ',alfaxz3_av(nclouds_xz)
- print*,'cp_useds1 ',(cp_useds1(icp),icp=1,ncp_tot)
- print*,'cp_useds2 ',(cp_useds2(icp),icp=1,ncp_tot)
+ print*,'- 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)
@@ -2402,7 +2790,7 @@
goto 91
endif
- if(DEBUG)then
+ if(DEBUG.EQ.1)then
print*,'---------------------- '
print*,'X-Z total clouds ',nclouds_xz
print*,' '
@@ -2454,10 +2842,9 @@
* -----------------------------------------------------------
* variables for track fitting
double precision AL_INI(5)
-c 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
@@ -2479,7 +2866,7 @@
enddo
enddo
ncp_ok=0
- do icp=1,ncp_tot !loop on couples
+ do icp=1,ncp_tot !loop over couples
* get info on
cpintersec(icp)=min(
$ cpcloud_yz(iyz,icp),
@@ -2488,6 +2875,7 @@
$ (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
+* cpintersec is >0 if yz and xz clouds contain the same image of couple icp
if(cpintersec(icp).ne.0)then
ncp_ok=ncp_ok+1
@@ -2520,16 +2908,18 @@
nplused=nplused+ hit_plane(ip)
enddo
-c if(nplused.lt.nplxz_min)goto 888 !next doublet
- if(nplused.lt.nplyz_min)goto 888 !next doublet
- if(ncp_ok.lt.ncpok)goto 888 !next cloud
- if(DEBUG)then
+ 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))
@@ -2575,7 +2965,8 @@
c$$$
c$$$ if(AL_INI(5).gt.defmax)goto 888 !next cloud
- if(DEBUG)then
+ 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))
@@ -2608,10 +2999,35 @@
hit_plane(6)=icp6
if(ncp_match(6).eq.0)hit_plane(6)=0 !-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)
@@ -2655,11 +3071,11 @@
ifail=0 !error flag in chi^2 computation
jstep=0 !number of minimization steps
iprint=0
-c if(DEBUG)iprint=1
- if(DEBUG)iprint=2
+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 *** '
$ //'(clouds_to_ctrack)'
@@ -2684,14 +3100,15 @@
* --------------------------
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
@@ -2699,7 +3116,7 @@
ntracks = ntracks + 1
- do ip=1,nplanes
+ do ip=1,nplanes !top to bottom
XV_STORE(ip,ntracks)=sngl(xv(ip))
YV_STORE(ip,ntracks)=sngl(yv(ip))
@@ -2716,6 +3133,7 @@
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))
@@ -2731,9 +3149,9 @@
SENSOR_STORE(nplanes-ip+1,ntracks)=0
LADDER_STORE(nplanes-ip+1,ntracks)=0
endif
- BX_STORE(nplanes-ip+1,ntracks)=0!I dont need it now
- BY_STORE(nplanes-ip+1,ntracks)=0!I dont need it now
- 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
@@ -2762,7 +3180,7 @@
return
endif
-c$$$ if(DEBUG)then
+c$$$ if(DEBUG.EQ.1)then
c$$$ print*,'****** TRACK CANDIDATES ***********'
c$$$ print*,'# R. chi2 RIG'
c$$$ do i=1,ntracks
@@ -2771,7 +3189,7 @@
c$$$ enddo
c$$$ print*,'***********************************'
c$$$ endif
- if(DEBUG)then
+ if(DEBUG.EQ.1)then
print*,'****** TRACK CANDIDATES *****************'
print*,'# R. chi2 RIG ndof'
do i=1,ntracks
@@ -2823,6 +3241,7 @@
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
@@ -2904,7 +3323,7 @@
LADDER_STORE(nplanes-ip+1,IBEST)=nldt
* --------------------------------------------------------------
- if(DEBUG)then
+ if(DEBUG.EQ.1)then
print*,
$ '------ Plane ',ip,' intersected on LADDER ',nldt
$ ,' SENSOR ',ist
@@ -2915,7 +3334,7 @@
* ===========================================
* STEP 1 >>>>>>> try to include a new couple
* ===========================================
-c if(DEBUG)print*,'>>>> try to include a new couple'
+c if(DEBUG.EQ.1)print*,'>>>> try to include a new couple'
distmin=1000000.
xmm = 0.
ymm = 0.
@@ -2947,7 +3366,7 @@
distance = distance_to(XP,YP)
c distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
id=id_cp(ip,icp,ist)
- if(DEBUG)print*,'( couple ',id
+ if(DEBUG.EQ.1)print*,'( couple ',id
$ ,' ) distance ',distance
if(distance.lt.distmin)then
xmm = xPAM
@@ -2979,7 +3398,7 @@
dedxtrk_y(nplanes-ip+1) = dedxmmy !<<<
* -----------------------------------
CP_STORE(nplanes-ip+1,ibest)=idm
- if(DEBUG)print*,'%%%% included couple ',idm
+ if(DEBUG.EQ.1)print*,'%%%% included couple ',idm
$ ,' (dist.= ',distmin,', cut ',clinc,' )'
goto 133 !next plane
endif
@@ -2987,7 +3406,7 @@
* STEP 2 >>>>>>> try to include a single cluster
* either from a couple or single
* ================================================
-c if(DEBUG)print*,'>>>> try to include a new cluster'
+c if(DEBUG.EQ.1)print*,'>>>> try to include a new cluster'
distmin=1000000.
xmm_A = 0. !---------------------------
ymm_A = 0. ! init variables that
@@ -3024,7 +3443,7 @@
$ )
distance = distance_to(XP,YP)
c distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
- if(DEBUG)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
@@ -3057,7 +3476,7 @@
$ )
distance = distance_to(XP,YP)
c distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
- if(DEBUG)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
@@ -3102,10 +3521,10 @@
distance = distance_to(XP,YP)
c distance = distance / RCHI2_STORE(ibest)!<<< MS !QUIQUI
- if(DEBUG)print*,'( cl-s ',icl
- $ ,' ) distance ',distance,'<',distmin,' ?'
+ if(DEBUG.EQ.1)print*,'( cl-s ',icl
+ $ ,' ) distance ',distance
if(distance.lt.distmin)then
- if(DEBUG)print*,'YES'
+c if(DEBUG.EQ.1)print*,'YES'
xmm_A = xPAM_A
ymm_A = yPAM_A
zmm_A = zPAM_A
@@ -3151,14 +3570,14 @@
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
+ if(DEBUG.EQ.1)print*,'%%%% included X-cl ',iclm
$ ,'( chi^2, ',RCHI2_STORE(ibest)
$ ,', dist.= ',distmin
$ ,', cut ',clinc,' )'
else
YGOOD(nplanes-ip+1)=1.
resy(nplanes-ip+1)=rymm
- if(DEBUG)print*,'%%%% included Y-cl ',iclm
+ if(DEBUG.EQ.1)print*,'%%%% included Y-cl ',iclm
$ ,'( chi^2, ',RCHI2_STORE(ibest)
$ ,', dist.= ', distmin
$ ,', cut ',clinc,' )'
@@ -3200,6 +3619,8 @@
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)
@@ -3208,10 +3629,10 @@
if(id.ne.0)then
iclx=clx(ip,icp_cp(id))
icly=cly(ip,icp_cp(id))
- cl_used(iclx)=ntrk !tag used clusters
- cl_used(icly)=ntrk !tag used clusters
+c$$$ cl_used(iclx)=ntrk !tag used clusters
+c$$$ cl_used(icly)=ntrk !tag used clusters
elseif(icl.ne.0)then
- cl_used(icl)=ntrk !tag used clusters
+c$$$ cl_used(icl)=ntrk !tag used clusters
endif
* -----------------------------
@@ -3230,7 +3651,7 @@
$ cly(ip,icp).eq.icl
$ )then
id=id_cp(ip,icp,1)
- if(DEBUG)then
+ if(DEBUG.EQ.1)then
print*,ip,' <<< cp ',id
$ ,' ( cl-x '
$ ,clx(ip,icp)
@@ -3490,12 +3911,12 @@
zv_nt(ip,ntr) = sngl(zv(ip))
axv_nt(ip,ntr) = sngl(axv(ip))
ayv_nt(ip,ntr) = sngl(ayv(ip))
-c l'avevo dimenticato!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
factor = sqrt(
$ tan( acos(-1.) * sngl(axv(ip)) /180. )**2 +
$ tan( acos(-1.) * sngl(ayv(ip)) /180. )**2 +
$ 1. )
-c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
dedx_x(ip,ntr) = sngl(dedxtrk_x(ip)/factor)
dedx_y(ip,ntr) = sngl(dedxtrk_y(ip)/factor)
@@ -3526,17 +3947,16 @@
cltrx(ip,ntr) = clx(nplanes-ip+1,icp_cp(id))
cltry(ip,ntr) = cly(nplanes-ip+1,icp_cp(id))
-c$$$ if(is_cp(id).ne.ssensor)
-c$$$ $ print*,'ERROR is sensor assignment (couple)'
-c$$$ $ ,is_cp(id),ssensor
-c$$$ if(LADDER(clx(nplanes-ip+1,icp_cp(id))).ne.sladder)
-c$$$ $ print*,'ERROR is ladder assignment (couple)'
-c$$$ $ ,LADDER(clx(nplanes-ip+1,icp_cp(id))),sladder
-
- nnnnx = npfastrips(clx(nplanes-ip+1,icp_cp(id)),PFA,angx)
- nnnny = npfastrips(cly(nplanes-ip+1,icp_cp(id)),PFA,angy)
- xbad(ip,ntr)= nbadstrips(nnnnx,clx(nplanes-ip+1,icp_cp(id)))
- ybad(ip,ntr)= nbadstrips(nnnny,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
+
+c$$$ nnnnx = npfastrips(clx(nplanes-ip+1,icp_cp(id)),PFA,angx)
+c$$$ nnnny = npfastrips(cly(nplanes-ip+1,icp_cp(id)),PFA,angy)
+c$$$ xbad(ip,ntr)= nbadstrips(nnnnx,clx(nplanes-ip+1,icp_cp(id)))
+c$$$ ybad(ip,ntr)= nbadstrips(nnnny,cly(nplanes-ip+1,icp_cp(id)))
+ xbad(ip,ntr)= nbadstrips(4,clx(nplanes-ip+1,icp_cp(id)))
+ ybad(ip,ntr)= nbadstrips(4,cly(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)
@@ -3544,25 +3964,35 @@
$ dedx_y(ip,ntr)=-dedx_y(ip,ntr)
elseif(icl.ne.0)then
-c >>> is a singlet
-c$$$ if(LADDER(icl).ne.sladder)
-c$$$ $ print*,'ERROR is ladder assignment (single)'
-c$$$ $ ,LADDER(icl),sladder
+
+ cl_used(icl) = 1 !tag used clusters
+
if(mod(VIEW(icl),2).eq.0)then
cltrx(ip,ntr)=icl
- nnnnn = npfastrips(icl,PFA,angx)
- xbad(ip,ntr) = nbadstrips(nnnnn,icl)
+c$$$ nnnnn = npfastrips(icl,PFA,angx)
+c$$$ xbad(ip,ntr) = nbadstrips(nnnnn,icl)
+ xbad(ip,ntr) = nbadstrips(4,icl)
+
if(nsatstrips(icl).gt.0)dedx_x(ip,ntr)=-dedx_x(ip,ntr)
elseif(mod(VIEW(icl),2).eq.1)then
cltry(ip,ntr)=icl
- nnnnn = npfastrips(icl,PFA,angy)
- ybad(ip,ntr) = nbadstrips(nnnnn,icl)
+c$$$ nnnnn = npfastrips(icl,PFA,angy)
+c$$$ ybad(ip,ntr) = nbadstrips(nnnnn,icl)
+ ybad(ip,ntr) = nbadstrips(4,icl)
if(nsatstrips(icl).gt.0)dedx_y(ip,ntr)=-dedx_y(ip,ntr)
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
+ endif
c$$$ print*,(xgood(i),i=1,6)
c$$$ print*,(ygood(i),i=1,6)
@@ -3593,12 +4023,19 @@
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
@@ -3638,6 +4075,30 @@
***** 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
+* --------------------------------------------------
+
+
+c$$$ print*,' cl ',icl,' --> ',cl_used(icl)
+c$$$
+c$$$ if( cl_used(icl).ne.0 )then
+c$$$ if(
+c$$$ $ mod(VIEW(icl),2).eq.0.and.
+c$$$ $ cltrx(ip,whichtrack(icl)).ne.icl )
+c$$$ $ print*,'**WARNING** cltrx(',ip,',',whichtrack(icl)
+c$$$ $ ,')=',cltrx(ip,whichtrack(icl)),'.ne.',icl
+c$$$ if(
+c$$$ $ mod(VIEW(icl),2).eq.1.and.
+c$$$ $ cltry(ip,whichtrack(icl)).ne.icl )
+c$$$ $ print*,'**WARNING** cltry(',ip,',',whichtrack(icl)
+c$$$ $ ,')=',cltry(ip,whichtrack(icl)),'.ne.',icl
+c$$$ endif
+
enddo
end