MODULE TPCtrack USE IOunits USE Threshold USE TPCdata USE TPCpaw USE MyMinuit USE Space PRIVATE ! eventVeto(iGroupV,iRow1,iRow2,iRowT,IVETO) ! findSeedTrack(iRow1,iRow2,iRowT,seedTrack,IERR) ! deMultiplex(seedTrack) PUBLIC :: eventVeto,findSeedTrack,deMultiplex ! drift time fits ! TimeFit(iRow1,iRow2,iRowT,time0,deltaT) staight line fit ! TimeLine(iRowT,time0,deltaT) linear regression ! TrackFitZ(iRowT,ZTrack) returns line, needs drift velocity etc. PUBLIC :: TimeLine,TimeFit,TrackFitZ ! spatial fits ! Track can be XYTrack_type: the radius is fit ! Line_type: straight track ! TrackFitWidth(iRowT,Track,Width) returns track and width ! TrackFitXY(iRowT,Z,Track) returns track, width is calculated from Z ! which can be ZTrack or T0 ! TrackFitX(iRowT,Z,Track,reso) returns resolution ! Z can be ZTrack or Width PUBLIC :: TrackFitWidth,TrackFitXY,TrackFitX ! trackWidth(ZTrack,y) function to calculate track width from ZTrack and Y PUBLIC :: trackWidth ! PRF(fitTrack,width,iPad) pad response function; fitTrack = Line ! ErrorOnAmplitude = AERR(Amplitude) PUBLIC :: PRF,AERR ! ==================================================================== PRIVATE :: GetStart,HitInRow,TrackHitInRow PRIVATE :: deltaT2Z,T2Z,trackWidthT,trackWidthZ PRIVATE :: ZCHI,TRKCHI,TRKCHIZ,ROWCHI,ROWCHItrunc,ROWCHIPRF PRIVATE :: TrackFitWidthRadius,TrackFitWidthLine PRIVATE :: TrackFitXYRadiusZ,TrackFitXYRadiusT,TrackFitXYLineZ,TrackFitXYLineT PRIVATE :: TrackFitXRadiusZ,TrackFitXRadiusW,TrackFitXLineZ,TrackFitXLineW PRIVATE :: QRAT,countHit INTERFACE TrackFitWidth MODULE PROCEDURE TrackFitWidthRadius,TrackFitWidthLine END INTERFACE INTERFACE TrackFitXY MODULE PROCEDURE TrackFitXYRadiusZ,TrackFitXYRadiusT, & TrackFitXYLineZ, TrackFitXYLineT END INTERFACE INTERFACE TrackFitX MODULE PROCEDURE TrackFitXRadiusZ,TrackFitXRadiusW, & TrackFitXLineZ, TrackFitXLineW END INTERFACE INTEGER,PRIVATE,SAVE,ALLOCATABLE,DIMENSION(:) :: iRowFit TYPE(Line_type),PRIVATE,SAVE :: ZTrackFit CONTAINS ! >>>> VETO >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> SUBROUTINE eventVeto(iGroupV,iRow1,iRow2,iRowT,IVETO) ! ************************************************************** ! *** iGroupV: array of groups, if one has a hit ! *** with amplitude > minSignalVeto IVETO = 1 ! *** iRow1,iRow2: arrays of seedRows, if no hit ! with amplitude > minSignalRow IVETO = 2 ! *** iRowT: array of tracking Rows, ! *** less than minHitTrack hits IVETO = 3 ! ************************************************************** INTEGER,DIMENSION(:),INTENT(IN) :: iGroupV,iRow1,iRow2,iRowT INTEGER,INTENT(OUT) :: IVETO INTEGER :: iG,iGroup,iR,nHit REAL :: Amplitude,T0 LOGICAL :: OK1,OK2 IVETO = 0 ! no hits in veto groups DO iG=1,SIZE(iGroupV) iGroup = iGroupV(iG) CALL GetSignalGroup(iGroup,Amplitude,T0) IF( Amplitude >= minSignalVeto )THEN IVETO = 1 RETURN ENDIF ENDDO ! hits in seed-rows OK1 = .FALSE. DO iR=1,SIZE(iRow1) IF( hitInRow(iRow1(iR),minSignalRow) ) THEN OK1 = .TRUE. EXIT ENDIF ENDDO OK2 = .FALSE. DO iR=1,SIZE(iRow2) IF( hitInRow(iRow2(iR),minSignalRow) )THEN OK2 = .TRUE. EXIT ENDIF ENDDO IF( .NOT.(OK1.AND.OK2) )THEN IVETO = 2 RETURN ENDIF ! hits in tracking-rows nHit = 0 DO iR=1,SIZE(iRowT) IF( hitInRow(iRowT(iR),minSignalRow) ) nHit = nHit + 1 ENDDO IF( nHit < minHitTrack )THEN IVETO = 3 RETURN ENDIF END SUBROUTINE eventVeto FUNCTION hitInRow(iRow,minAmpl) RESULT(OK) INTEGER,INTENT(IN) :: iRow REAL,INTENT(IN) :: minAmpl LOGICAL :: OK INTEGER :: iP,iPad REAL :: Amplitude,T0 OK = .FALSE. DO iP=1,nPadsinRow(iRow) iPad = iPadofRow(iP,iRow) CALL GetSignal(iPad,Amplitude,T0) IF( Amplitude >= minAmpl )THEN OK = .TRUE. RETURN ENDIF ENDDO END FUNCTION hitInRow FUNCTION trackHitInRow(iRow,myLine,minAmpl) RESULT(OK) INTEGER,INTENT(IN) :: iRow TYPE(Line_type),INTENT(IN) :: myLine REAL,INTENT(IN) :: minAmpl LOGICAL :: OK INTEGER :: iP,iPad REAL :: Amplitude,T0,dist OK = .FALSE. DO iP=1,nPadsinRow(iRow) iPad = iPadofRow(iP,iRow) CALL GetSignal(iPad,Amplitude,T0) IF( Amplitude >= minAmpl )THEN dist = Distance(PadLocation(iPad),myLine) ! IF( dist <= 0.0 )THEN : track must cross pad ! IF( dist <= x.x )THEN : track must be closer than x.x mm to pad IF( dist <= 0.0 )THEN OK = .TRUE. RETURN ENDIF ENDIF ENDDO END FUNCTION trackHitInRow ! <<<< VETO <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< ! >>>> TRACKING >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> SUBROUTINE findSeedTrack(iRow1,iRow2,iRowT,seedTrack,IERR) ! ************************************************************** ! *** iRow1,iRow2: arrays of seedRows ! *** iRowT: array of tracking rows ! *** find largest signal in iRow1 and iRow2, connect with line ! *** seedTrack: line with most hits in iRowT ! *** IERR: number of seedTracks-1 => 0 is OK ! ************************************************************** USE Space INTEGER,DIMENSION(:),INTENT(IN) :: iRow1,iRow2,iRowT TYPE(Line_type),INTENT(OUT) :: seedTrack TYPE(Line_type) :: myTrack INTEGER,INTENT(OUT) :: IERR TYPE(Coordinate_type) :: Point1,Point2 INTEGER :: minHit,nTrack,nHitRow,iGroup,iGroup1,iGroup2 INTEGER :: iR,iP1,iPad1,iP2,iPad2 REAL :: maxAmpl,ampl,T0 LOGICAL :: OK1,OK2 IERR = -3 nTrack = 0 minHit = minHitTrack ! get group with maximum amplitude in both sets of seed-rows maxAmpl = 0.0 iGroup1 = -1 DO iR=1,SIZE(iRow1) CALL GmaxSignalRow(iRow1(iR),iGroup,ampl,T0) IF( ampl > maxAmpl )THEN maxAmpl = ampl iGroup1 = iGroup ENDIF ENDDO maxAmpl = 0.0 iGroup2 = -1 DO iR=1,SIZE(iRow2) CALL GmaxSignalRow(iRow2(iR),iGroup,ampl,T0) IF( ampl > maxAmpl )THEN maxAmpl = ampl iGroup2 = iGroup ENDIF ENDDO IF( iGroup1<1 .OR. iGroup2<1 )THEN IERR = -2 RETURN ENDIF ! loop over pads of these groups within the rows ROW1: DO iP1 = 1,nPadsinGroup(iGroup1) iPad1 = iPadofGroup(iP1,iGroup1) ! is this pad in one of the rows? OK1 = .FALSE. DO iR=1,SIZE(iRow1) IF( PadRow(iPad1) == iRow1(iR) ) THEN OK1 = .TRUE. EXIT ENDIF ENDDO IF( OK1 )THEN CALL GetSignal(iPad1,ampl,T0) IF( ampl < minSignalRow )THEN IERR = -1 RETURN ENDIF CALL getStart(iPad1,iRow1,Point1) ROW2: DO iP2 = 1,nPadsinGroup(iGroup2) iPad2 = iPadofGroup(iP2,iGroup2) ! is this pad in one of the rows? OK2 = .FALSE. DO iR=1,SIZE(iRow2) IF( PadRow(iPad2) == iRow2(iR) ) THEN OK2 = .TRUE. EXIT ENDIF ENDDO IF( OK2 )THEN CALL GetSignal(iPad2,ampl,T0) IF( ampl < minSignalRow )THEN IERR = -1 RETURN ENDIF CALL getStart(iPad2,iRow2,Point2) myTrack = Line(Point1,Point2) ! count hits in rows nHitRow = 0 DO iR=1,SIZE(iRowT) IF( trackHitInRow(iRowT(iR),myTrack,minSignalRow) ) & nHitRow = nHitRow+1 ENDDO ! do we have a new seedtrack? IF( nHitRow == minHit )THEN seedTrack = myTrack nTrack = nTrack+1 ELSEIF( nHitRow > minHit )THEN minHit = nHitRow seedTrack = myTrack ! nTrack = nTrack+1 nTrack = 1 ENDIF ENDIF ENDDO ROW2 ENDIF ENDDO ROW1 ! returns IERR=0 if exactly one track is found IERR = nTrack - 1 END SUBROUTINE findSeedTrack SUBROUTINE getStart(iPad,iRow,Point) ! gets a start value for the track position ! uses this pad and it's neighbours within rows iRow ! and calculates the average position INTEGER,INTENT(IN) :: iPad INTEGER,DIMENSION(:),INTENT(IN) :: iRow TYPE(Coordinate_type),INTENT(OUT) :: Point INTEGER :: i,iR,j,iP REAL :: Amplitude,T0 TYPE(Coordinate_type) :: Pkt REAL(KIND=spPrec) :: sumAmpl,dist CALL GetSignal(iPad,Amplitude,T0) sumAmpl = Amplitude Pkt = getCentre(PadLocation(iPad)) Point = Pkt*Amplitude ! loop over adjacent pads in seedrow DO i = 1,SIZE(iRow) iR = iRow(i) DO j = 1,nPadsinRow(iR) iP = iPadofRow(j,iR) IF( iP /= iPad )THEN dist = Distance( PadLocation(iPad) , PadLocation(iP) ) IF( dist < minDistPad )THEN CALL GetSignal(iP,Amplitude,T0) sumAmpl = sumAmpl + Amplitude Pkt = getCentre(PadLocation(iPad)) Point = Point + Pkt*Amplitude ENDIF ENDIF ENDDO ENDDO IF( sumAmpl>0.0 )THEN sumAmpl = 1.0 / sumAmpl Point = Point * sumAmpl ELSE Point = Coordinate(-1000.0,-1000.0) ENDIF END SUBROUTINE getStart SUBROUTINE deMultiplex(seedTrack) ! ************************************************************** ! *** set all amplitudes to -1 ! *** apart from one pad in each group closest to seedTrack ! ************************************************************** TYPE(Line_type),INTENT(IN) :: seedTrack INTEGER :: iGroup,iPadNear,iP,iPad REAL(KIND=spPrec) :: minDist,dist DO iGroup=1,nGroups minDist = 1.0E6 iPadNear = -1 DO iP=1,nPadsinGroup(iGroup) iPad = iPadofGroup(iP,iGroup) dist = Distance(PadLocation(iPad),seedTrack) IF( dist < minDist )THEN minDist = dist iPadNear = iPad ENDIF ENDDO DO iP=1,nPadsinGroup(iGroup) iPad = iPadofGroup(iP,iGroup) IF( iPad /= iPadNear ) CALL VoidSignal(iPad) ENDDO ENDDO END SUBROUTINE deMultiplex FUNCTION T2Z(T0) RESULT(DD) ! get driftdistance from timebin T0 ! T0: drift time in time bins ! DD: drift distance in mm REAL,INTENT(IN) :: T0 REAL :: DD DD = T0 IF( tBin2ns<=0.0 .OR. vDrift<=0.0 )THEN PRINT*,"ERROR T2Z> parameters not given" RETURN ENDIF DD = deltaT2Z(T0-tBin0) END FUNCTION T2Z FUNCTION deltaT2Z(T0) RESULT(DD) ! get driftdistance from timebin difference T0 ! T0: drift time in time bins ! DD: drift distance in mm REAL,INTENT(IN) :: T0 REAL :: DD DD = T0 IF( tBin2ns<=0.0 .OR. vDrift<=0.0 )THEN PRINT*,"ERROR deltaT2Z> parameters not given" RETURN ENDIF ! vDrift is given in cm/ns; convert to mm DD = 10.0 * T0*tBin2ns * vDrift END FUNCTION deltaT2Z FUNCTION trackWidthT(T0) RESULT(tWidth) ! get driftdistance from timebin T0 REAL,INTENT(IN) :: T0 REAL :: DD REAL :: tWidth DD = T2Z(T0) tWidth = trackWidthZ(DD) END FUNCTION trackWidthT FUNCTION trackWidthZ(DD) RESULT(tWidth) ! get track width from trackWidthPar ! driftdistance DD in mm, trackwidth in mm REAL,INTENT(IN) :: DD REAL :: D,tWidth tWidth = -1.0 IF( size(trackWidthPar)<1 )THEN PRINT*,"ERROR trackWidth> parameters not given" RETURN ENDIF IF( DD>0.0 )THEN ! trackWidthPar for driftDistance in cm, convert D = 0.1 * DD tWidth = POLY(D,trackWidthPar) tWidth = sqrt( tWidth ) ELSE tWidth = sqrt( trackWidthPar(1) ) ENDIF END FUNCTION trackWidthZ FUNCTION trackWidth(ZTrack,y) RESULT(tWidth) ! get driftdistance from Ztrack TYPE(Line_type),INTENT(IN) :: ZTrack REAL,INTENT(IN) :: Y REAL :: DD,tWidth DD = getLineX(ZTrack,Y) tWidth = trackWidthZ(DD) END FUNCTION trackWidth ! <<<< TRACKING <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< ! >>>> MINUIT >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> SUBROUTINE TrackFitZ(iRowT,ZTrack) INTEGER,DIMENSION(:),INTENT(IN) :: iRowT TYPE(Line_type),INTENT(OUT) :: ZTrack REAL :: time0, deltaT, deltaY=100.0 ! deltaY is any number REAL :: Z,dZ TYPE(Coordinate_type) :: A,B CALL TimeLine(iRowT,time0, deltaT) IF( time0<-999.0 .and. deltaT<-999.0 )THEN ZTrack = Line(0.0,0.0) ELSE Z = T2Z(time0) dZ = Z + deltaT2Z(deltaT)*deltaY A = Coordinate(Z,0.0) B = Coordinate(dZ,deltaY) ZTrack = Line(A,B) ENDIF END SUBROUTINE TrackFitZ SUBROUTINE TimeLine(iRowT,time0,deltaT) ! TimeFit via LinearRegression ! Error: time0 = -1000.0, deltaT = -1000.0 INTEGER,DIMENSION(:),INTENT(IN) :: iRowT REAL,INTENT(OUT) :: time0, deltaT REAL,DIMENSION(:),ALLOCATABLE :: X,Y,SIGSQ REAL :: AmpRow, Trow INTEGER :: I, iRow, iPad, IERR, iGmaxAmpl ALLOCATE(X(size(iRowT)),STAT=IERR) ALLOCATE(Y(size(iRowT)),STAT=IERR) ALLOCATE(SIGSQ(size(iRowT)),STAT=IERR) time0 = 0.0 deltaT = 0.0 X = 0.0 Y = 0.0 SIGSQ = 0.0 DO I=1,size(iRowT) iRow = iRowT(I) ! CALL GetSignalRow(iRow, AmpRow, Trow) CALL GmaxSignalRow(iRow,iGmaxAmpl, AmpRow, Trow) IF( AmpRow < minSignal ) CYCLE iPad = iPadofRow(1,iRow) X(I) = Length(getCentre(PadLocation(iPad))) Y(I) = Trow SIGSQ(I) = 1.0 ENDDO CALL LinearRegression(X,Y,SIGSQ,time0,deltaT,IERR) IF( IERR /= 0 )THEN time0 = -1000.0 deltaT = -1000.0 ENDIF DEALLOCATE(X) DEALLOCATE(Y) DEALLOCATE(SIGSQ) END SUBROUTINE TimeLine ! track fit y-z 'theta' ******************************************* SUBROUTINE TimeFit(iRow1,iRow2,iRowT,time0,deltaT) INTEGER,DIMENSION(:),INTENT(IN) :: iRow1,iRow2,iRowT REAL,INTENT(OUT) :: time0,deltaT CHARACTER(LEN=50) :: CHTITL INTEGER,PARAMETER :: NFITP=2 CHARACTER(LEN=10),DIMENSION(NFITP) :: PCH REAL(KIND=minuitPrec),DIMENSION(NFITP) :: PVAL,PERR,BND1,BND2 REAL(KIND=minuitPrec),DIMENSION(NFITP,NFITP) :: PMAT REAL(KIND=minuitPrec) :: CHI2 REAL :: Amplitude,T1,T2,Y1,Y2 INTEGER :: IERR,IPRT,iPad ! get the start values from the seed rows CALL GetSignalRow(iRow1(1),Amplitude,T1) CALL GetSignalRow(iRow2(1),Amplitude,T2) iPad = iPadofRow(1,iRow1(1)) Y1 = Length(getCentre(PadLocation(iPad))) iPad = iPadofRow(1,iRow2(1)) Y2 = Length(getCentre(PadLocation(iPad))) PCH(1) = "T0 " PCH(2) = "TSLOPE " PVAL(2) = ( T1-T2 )/( Y1-Y2 ) PVAL(1) = T1 - PVAL(2)*Y1 PERR(1) = 2.0 PERR(2) = 0.3 ! normaly no bounds BND1 = 0.0 BND2 = 0.0 ! define the rows for the fit ALLOCATE( iRowFit(size(iRowT)),STAT=IERR ) iRowFit = iRowT CHTITL = "+++++ TPC Z fit +++++" IPRT = 0 CALL MINUIT(ZCHI,CHTITL,NFITP,PCH,PVAL,PERR,PMAT, & BND1,BND2,CHI2,IPRT,IERR) DEALLOCATE(iRowFit) IF( IERR == 0 .OR. IERR == 1 )THEN time0 = PVAL(1) deltaT = PVAL(2) ! print"(3F12.6)",PERR,CHI2 ELSE time0 = -1000.0 deltaT = -1000.0 ENDIF END SUBROUTINE TimeFit SUBROUTINE ZCHI(NPAR,GRAD,FVAL,XVAL,IFLAG) ! subroutine to be called by MINUIT - calculates the chisq INTEGER,INTENT(IN) :: IFLAG , NPAR REAL(KIND=minuitPrec),DIMENSION(2),INTENT(IN) :: XVAL REAL(KIND=minuitPrec),DIMENSION(2),INTENT(OUT) :: GRAD REAL(KIND=minuitPrec),INTENT(OUT) :: FVAL REAL(KIND=minuitPrec) :: CHISQ INTEGER :: I,iRow,iPad,myPar REAL :: Trow,AmpRow REAL(KIND=minuitPrec) :: T,Y,Amplitude,Tline REAL(KIND=minuitPrec),PARAMETER :: TERR=20.0 ! dummy to avoid compiler warnings IF( IFLAG==9 ) myPar=NPAR GRAD = 0.0_minuitPrec CHISQ = 0.0_minuitPrec DO I=1,size(iRowFit) iRow = iRowFit(I) CALL GetSignalRow(iRow,AmpRow,Trow) IF( Trow<0.0 )CYCLE T = Trow Amplitude = AmpRow iPad = iPadofRow(1,iRow) Y = Length(getCentre(PadLocation(iPad))) Tline = XVAL(1) + XVAL(2)*Y CHISQ = CHISQ + ((T-Tline)/TERR)**2 ENDDO FVAL = CHISQ END SUBROUTINE ZCHI ! TrackFit *************************************************** SUBROUTINE TrackFitWidthRadius(iRowT,Track,Width) ! fit 4 parameters: X0, phi, invRad, width INTEGER,DIMENSION(:),INTENT(IN) :: iRowT TYPE(XYTrack_type),INTENT(INOUT) :: Track REAL(KIND=spPrec),INTENT(OUT) :: Width CHARACTER(LEN=50) :: CHTITL INTEGER,PARAMETER :: NFITP=4 CHARACTER(LEN=10),DIMENSION(NFITP) :: PCH REAL(KIND=minuitPrec),DIMENSION(NFITP) :: PVAL,PERR,BND1,BND2 REAL(KIND=minuitPrec),DIMENSION(NFITP,NFITP) :: PMAT REAL(KIND=minuitPrec) :: CHI2 REAL(KIND=spPrec) :: X0,Phi,invRad INTEGER :: IERR,IPRT CALL getXYTrack(Track,X0,Phi,invRad) PCH(1) = "X0 " PCH(2) = "PHI " PCH(3) = "INVRAD " PCH(4) = "SIGMA " PVAL(1) = X0 ! to avoid phi running away, fit tan(phi) PVAL(2) = tan(Phi) PVAL(3) = invRad PVAL(4) = 1.0 PERR(1) = 0.2 PERR(2) = 0.1 PERR(3) = 0.01 PERR(4) = 0.3 ! normaly no bounds BND1 = 0.0 BND2 = 0.0 ! define the ZTrack ZTrackFit = Line(0.0,0.0) ! define the rows for the fit ALLOCATE( iRowFit(size(iRowT)),STAT=IERR ) iRowFit = iRowT CHTITL = "+++++ TPC track fit, with radius and width +++++" IPRT = 0 CALL MINUIT(TRKCHI,CHTITL,NFITP,PCH,PVAL,PERR,PMAT, & BND1,BND1,CHI2,IPRT,IERR) DEALLOCATE(iRowFit) IF( IERR == 0 .OR. IERR == 1 )THEN X0 = PVAL(1) Phi = atan(PVAL(2)) invRad = PVAL(3) Width = PVAL(4) ELSE WRITE(UNIT=IOErr,FMT=*) "WARNING,TRKFIT > track fit failed:" & ,IERR," Event:",EventNumber ! keep the seed track parameters Width = 1000.0 ENDIF Track = XYTrack(Line(X0,Phi),invRad) END SUBROUTINE TrackFitWidthRadius SUBROUTINE TrackFitWidthLine(iRowT,Track,Width) ! fit 3 parameters: X0, phi, width INTEGER,DIMENSION(:),INTENT(IN) :: iRowT TYPE(Line_type),INTENT(INOUT) :: Track REAL(KIND=spPrec),INTENT(OUT) :: Width CHARACTER(LEN=50) :: CHTITL INTEGER,PARAMETER :: NFITP=4 CHARACTER(LEN=10),DIMENSION(NFITP) :: PCH REAL(KIND=minuitPrec),DIMENSION(NFITP) :: PVAL,PERR,BND1,BND2 REAL(KIND=minuitPrec),DIMENSION(NFITP,NFITP) :: PMAT REAL(KIND=minuitPrec) :: CHI2 REAL(KIND=spPrec) :: X0,Phi INTEGER :: IERR,IPRT CALL getLine(Track,X0,Phi) PCH(1) = "X0 " PCH(2) = "PHI " PCH(3) = "INVRAD " PCH(4) = "SIGMA " PVAL(1) = X0 PVAL(2) = tan(Phi) PVAL(3) = 0.0 PVAL(4) = 1.0 PERR(1) = 0.2 PERR(2) = 0.1 PERR(3) =-1.0 PERR(4) = 0.3 ! normaly no bounds BND1 = 0.0 BND2 = 0.0 ! define the ZTrack ZTrackFit = Line(0.0,0.0) ! define the rows for the fit ALLOCATE( iRowFit(size(iRowT)),STAT=IERR ) iRowFit = iRowT CHTITL = "+++++ TPC track fit , variable width +++++" IPRT = 0 CALL MINUIT(TRKCHI,CHTITL,NFITP,PCH,PVAL,PERR,PMAT, & BND1,BND1,CHI2,IPRT,IERR) DEALLOCATE(iRowFit) IF( IERR == 0 .OR. IERR == 1 )THEN X0 = PVAL(1) Phi = atan(PVAL(2)) Width = PVAL(4) ELSE WRITE(UNIT=IOErr,FMT=*) "WARNING,TRKFIT > track fit failed:" & ,IERR," Event:",EventNumber ! keep the seed track parameters Width = 1000.0 ENDIF Track = Line(X0,Phi) END SUBROUTINE TrackFitWidthLine ! TrackFitXY ****************************************************** ! track width is taken from input SUBROUTINE TrackFitXYRadiusZ(iRowT,ZTrack,Track) ! fit 3 parameters: X0, phi, invRad ! width is calculated from ZTrack INTEGER,DIMENSION(:),INTENT(IN) :: iRowT TYPE(Line_type),INTENT(IN) :: ZTrack TYPE(XYTrack_type),INTENT(INOUT) :: Track CHARACTER(LEN=50) :: CHTITL INTEGER,PARAMETER :: NFITP=3 CHARACTER(LEN=10),DIMENSION(NFITP) :: PCH REAL(KIND=minuitPrec),DIMENSION(NFITP) :: PVAL,PERR,BND1,BND2 REAL(KIND=minuitPrec),DIMENSION(NFITP,NFITP) :: PMAT REAL(KIND=minuitPrec) :: CHI2 REAL(KIND=spPrec) :: X0,Phi,invRad INTEGER :: IERR,IPRT CALL getXYTrack(Track,X0,Phi,invRad) PCH(1) = "X0 " PCH(2) = "PHI " PCH(3) = "INVRAD " PVAL(1) = X0 PVAL(2) = tan(Phi) PVAL(3) = invRad PERR(1) = 0.2 PERR(2) = 0.1 PERR(3) = 0.01 ! normaly no bounds BND1 = 0.0 BND2 = 0.0 ! define the ZTrack ZTrackFit = ZTrack ! define the rows for the fit ALLOCATE( iRowFit(size(iRowT)),STAT=IERR ) iRowFit = iRowT CHTITL = "+++++ TPC track fit, with radius +++++" IPRT = 0 CALL MINUIT(TRKCHIZ,CHTITL,NFITP,PCH,PVAL,PERR,PMAT, & BND1,BND1,CHI2,IPRT,IERR) DEALLOCATE(iRowFit) IF( IERR == 0 .OR. IERR == 1 )THEN X0 = PVAL(1) Phi = atan(PVAL(2)) invRad = PVAL(3) ELSE WRITE(UNIT=IOErr,FMT=*) "WARNING,TRKFIT > track fit failed:" & ,IERR," Event:",EventNumber ! keep the seed track parameters ENDIF Track = XYTrack(Line(X0,Phi),invRad) END SUBROUTINE TrackFitXYRadiusZ SUBROUTINE TrackFitXYRadiusT(iRowT,T0,Track) ! fit 3 parameters: X0, phi, invRad ! width is calculated from input INTEGER,DIMENSION(:),INTENT(IN) :: iRowT REAL,INTENT(IN) :: T0 TYPE(XYTrack_type),INTENT(INOUT) :: Track CHARACTER(LEN=50) :: CHTITL INTEGER,PARAMETER :: NFITP=4 CHARACTER(LEN=10),DIMENSION(NFITP) :: PCH REAL(KIND=minuitPrec),DIMENSION(NFITP) :: PVAL,PERR,BND1,BND2 REAL(KIND=minuitPrec),DIMENSION(NFITP,NFITP) :: PMAT REAL(KIND=minuitPrec) :: CHI2 REAL(KIND=spPrec) :: X0,Phi,invRad INTEGER :: IERR,IPRT CALL getXYTrack(Track,X0,Phi,invRad) PCH(1) = "X0 " PCH(2) = "PHI " PCH(3) = "INVRAD " PCH(4) = "SIGMA " PVAL(1) = X0 PVAL(2) = tan(Phi) PVAL(3) = invRad PVAL(4) = trackWidthT(T0) ! if trackwidth < 0 there is something wrong IF( PVAL(4) < 0.0 ) STOP PERR(1) = 0.2 PERR(2) = 0.1 PERR(3) = 0.01 PERR(4) = -1.0 ! normaly no bounds BND1 = 0.0 BND2 = 0.0 ! define the ZTrack ZTrackFit = Line(0.0,0.0) ! define the rows for the fit ALLOCATE( iRowFit(size(iRowT)),STAT=IERR ) iRowFit = iRowT CHTITL = "+++++ TPC track fit, with radius +++++" IPRT = 0 CALL MINUIT(TRKCHI,CHTITL,NFITP,PCH,PVAL,PERR,PMAT, & BND1,BND1,CHI2,IPRT,IERR) DEALLOCATE(iRowFit) IF( IERR == 0 .OR. IERR == 1 )THEN X0 = PVAL(1) Phi = atan(PVAL(2)) invRad = PVAL(3) ELSE WRITE(UNIT=IOErr,FMT=*) "WARNING,TRKFIT > track fit failed:" & ,IERR," Event:",EventNumber ! keep the seed track parameters ENDIF Track = XYTrack(Line(X0,Phi),invRad) END SUBROUTINE TrackFitXYRadiusT SUBROUTINE TrackFitXYLineZ(iRowT,ZTrack,Track) ! fit 2 parameters: X0, phi ! width is calculated from ZTrack INTEGER,DIMENSION(:),INTENT(IN) :: iRowT TYPE(Line_type),INTENT(IN) :: ZTrack TYPE(Line_type),INTENT(INOUT) :: Track CHARACTER(LEN=50) :: CHTITL INTEGER,PARAMETER :: NFITP=3 CHARACTER(LEN=10),DIMENSION(NFITP) :: PCH REAL(KIND=minuitPrec),DIMENSION(NFITP) :: PVAL,PERR,BND1,BND2 REAL(KIND=minuitPrec),DIMENSION(NFITP,NFITP) :: PMAT REAL(KIND=minuitPrec) :: CHI2 REAL(KIND=spPrec) :: X0,Phi INTEGER :: IERR,IPRT CALL getLine(Track,X0,Phi) PCH(1) = "X0 " PCH(2) = "PHI " PCH(3) = "INVRAD " PVAL(1) = X0 PVAL(2) = tan(Phi) PVAL(3) = 0.0 PERR(1) = 0.2 PERR(2) = 0.1 PERR(3) = -1.0 ! normaly no bounds BND1 = 0.0 BND2 = 0.0 ! define the ZTrack ZTrackFit = ZTrack ! define the rows for the fit ALLOCATE( iRowFit(size(iRowT)),STAT=IERR ) iRowFit = iRowT CHTITL = "+++++ TPC track fit, invRad = 0 +++++" IPRT = 0 CALL MINUIT(TRKCHIZ,CHTITL,NFITP,PCH,PVAL,PERR,PMAT, & BND1,BND1,CHI2,IPRT,IERR) DEALLOCATE(iRowFit) IF( IERR == 0 .OR. IERR == 1 )THEN X0 = PVAL(1) Phi = atan(PVAL(2)) ELSE WRITE(UNIT=IOErr,FMT=*) "WARNING,TRKFIT > track fit failed:" & ,IERR," Event:",EventNumber ! keep the seed track parameters ENDIF Track = Line(X0,Phi) END SUBROUTINE TrackFitXYLineZ SUBROUTINE TrackFitXYLineT(iRowT,T0,Track) ! fit 2 parameters: X0, phi ! width is taken from input INTEGER,DIMENSION(:),INTENT(IN) :: iRowT REAL,INTENT(IN) :: T0 TYPE(Line_type),INTENT(INOUT) :: Track CHARACTER(LEN=50) :: CHTITL INTEGER,PARAMETER :: NFITP=4 CHARACTER(LEN=10),DIMENSION(NFITP) :: PCH REAL(KIND=minuitPrec),DIMENSION(NFITP) :: PVAL,PERR,BND1,BND2 REAL(KIND=minuitPrec),DIMENSION(NFITP,NFITP) :: PMAT REAL(KIND=minuitPrec) :: CHI2 REAL(KIND=spPrec) :: X0,Phi INTEGER :: IERR,IPRT CALL getLine(Track,X0,Phi) PCH(1) = "X0 " PCH(2) = "PHI " PCH(3) = "INVRAD " PCH(4) = "SIGMA " PVAL(1) = X0 PVAL(2) = tan(Phi) PVAL(3) = 0.0 PVAL(4) = trackWidthT(T0) ! if trackwidth < 0 there is something wrong IF( PVAL(4) < 0.0 ) STOP PERR(1) = 0.2 PERR(2) = 0.1 PERR(3) = -1.0 PERR(4) = -1.0 ! normaly no bounds BND1 = 0.0 BND2 = 0.0 ! define the ZTrack ZTrackFit = Line(0.0,0.0) ! define the rows for the fit ALLOCATE( iRowFit(size(iRowT)),STAT=IERR ) iRowFit = iRowT CHTITL = "+++++ TPC track fit, width fixed +++++" IPRT = 0 CALL MINUIT(TRKCHI,CHTITL,NFITP,PCH,PVAL,PERR,PMAT, & BND1,BND1,CHI2,IPRT,IERR) DEALLOCATE(iRowFit) IF( IERR == 0 .OR. IERR == 1 )THEN X0 = PVAL(1) Phi = atan(PVAL(2)) ELSE WRITE(UNIT=IOErr,FMT=*) "WARNING,TRKFIT > track fit failed:" & ,IERR," Event:",EventNumber ! keep the seed track parameters ENDIF Track = Line(X0,Phi) END SUBROUTINE TrackFitXYLineT ! resolution fit, 1-parameter: x ************************************* SUBROUTINE TrackFitXRadiusZ(iRowT,ZTrack,Track,reso) INTEGER,DIMENSION(:),INTENT(IN) :: iRowT TYPE(Line_type),INTENT(IN) :: ZTrack TYPE(XYTrack_type),INTENT(IN) :: Track REAL,INTENT(OUT) :: reso CHARACTER(LEN=50) :: CHTITL INTEGER,PARAMETER :: NFITP=3 CHARACTER(LEN=10),DIMENSION(NFITP) :: PCH REAL(KIND=minuitPrec),DIMENSION(NFITP) :: PVAL,PERR,BND1,BND2 REAL(KIND=minuitPrec),DIMENSION(NFITP,NFITP) :: PMAT REAL(KIND=minuitPrec) :: CHI2 REAL(KIND=spPrec) :: X0,Phi,invRad,Xrow,Yrow INTEGER :: IERR,IPRT,nHit CALL getXYTrack(Track,X0,Phi,invRad) Yrow = Length(getCentre(PadLocation(iPadofRow(1,iRowT(1))))) CALL countHit(iRowT,Line(Track,Yrow),nHit,Xrow) IF( nHit == 0 )THEN reso = -1000.0 RETURN ELSEIF( nHit == 1 )THEN ! fit will fail, set position to center of pad reso = Xrow - X0 RETURN ENDIF PCH(1) = "X0 " PCH(2) = "PHI " PCH(3) = "INVRAD " PVAL(1) = X0 PVAL(2) = tan(Phi) PVAL(3) = invRad PERR(1) = 0.1 PERR(2) = -0.03 PERR(3) = -0.03 ! normaly no bounds BND1 = 0.0 BND2 = 0.0 ! define the ZTrack ZTrackFit = ZTrack ! define the rows for the fit ALLOCATE( iRowFit(size(iRowT)),STAT=IERR ) iRowFit = iRowT WRITE(UNIT=CHTITL,FMT=*) "+++++ TPC resolution fit, rows:",iRowFit IPRT = 0 CALL MINUIT(TRKCHIZ,CHTITL,NFITP,PCH,PVAL,PERR,PMAT, & BND1,BND2,CHI2,IPRT,IERR) DEALLOCATE(iRowFit) IF( IERR == 0 .OR. IERR == 1 )THEN reso = PVAL(1) - X0 ELSE WRITE(UNIT=IOErr,FMT=*) "WARNING,TRKFIT > resolution fit failed:" & ,IERR," Event:",EventNumber reso = 1000.0 ENDIF END SUBROUTINE TrackFitXRadiusZ SUBROUTINE TrackFitXRadiusW(iRowT,Width,Track,reso) INTEGER,DIMENSION(:),INTENT(IN) :: iRowT REAL(KIND=spPrec),INTENT(IN) :: Width TYPE(XYTrack_type),INTENT(IN) :: Track REAL,INTENT(OUT) :: reso CHARACTER(LEN=50) :: CHTITL INTEGER,PARAMETER :: NFITP=4 CHARACTER(LEN=10),DIMENSION(NFITP) :: PCH REAL(KIND=minuitPrec),DIMENSION(NFITP) :: PVAL,PERR,BND1,BND2 REAL(KIND=minuitPrec),DIMENSION(NFITP,NFITP) :: PMAT REAL(KIND=minuitPrec) :: CHI2 REAL(KIND=spPrec) :: X0,Phi,invRad,Xrow,Yrow INTEGER :: IERR,IPRT,nHit CALL getXYTrack(Track,X0,Phi,invRad) Yrow = Length(getCentre(PadLocation(iPadofRow(1,iRowT(1))))) CALL countHit(iRowT,Line(Track,Yrow),nHit,Xrow) IF( nHit == 0 )THEN reso = -1000.0 RETURN ELSEIF( nHit == 1 )THEN ! fit will fail, set position to center of pad reso = Xrow - X0 RETURN ENDIF ! is Trackin good? IF( width > 999.0 )THEN reso = -1000.0 RETURN ENDIF PCH(1) = "X0 " PCH(2) = "PHI " PCH(3) = "INVRAD " PCH(4) = "SIGMA " PVAL(1) = X0 PVAL(2) = tan(Phi) PVAL(3) = invRad PVAL(4) = Width PERR(1) = 0.1 PERR(2) = -0.03 PERR(3) = -1.0 PERR(4) = -0.3 ! normaly no bounds BND1 = 0.0 BND2 = 0.0 ! define the ZTrack ZTrackFit = Line(0.0,0.0) ! define the rows for the fit ALLOCATE( iRowFit(size(iRowT)),STAT=IERR ) iRowFit = iRowT WRITE(UNIT=CHTITL,FMT=*) "+++++ TPC resolution fit, rows:",iRowFit IPRT = 0 CALL MINUIT(TRKCHI,CHTITL,NFITP,PCH,PVAL,PERR,PMAT, & BND1,BND2,CHI2,IPRT,IERR) DEALLOCATE(iRowFit) IF( IERR == 0 .OR. IERR == 1 )THEN reso = PVAL(1) - X0 ELSE WRITE(UNIT=IOErr,FMT=*) "WARNING,TRKFIT > resolution fit failed:" & ,IERR," Event:",EventNumber reso = 1000.0 ENDIF END SUBROUTINE TrackFitXRadiusW SUBROUTINE TrackFitXLineZ(iRowT,ZTrack,Track,reso) INTEGER,DIMENSION(:),INTENT(IN) :: iRowT TYPE(Line_type),INTENT(IN) :: ZTrack TYPE(Line_type),INTENT(IN) :: Track REAL,INTENT(OUT) :: reso CHARACTER(LEN=50) :: CHTITL INTEGER,PARAMETER :: NFITP=3 CHARACTER(LEN=10),DIMENSION(NFITP) :: PCH REAL(KIND=minuitPrec),DIMENSION(NFITP) :: PVAL,PERR,BND1,BND2 REAL(KIND=minuitPrec),DIMENSION(NFITP,NFITP) :: PMAT REAL(KIND=minuitPrec) :: CHI2 REAL(KIND=spPrec) :: X0,Phi,Xrow INTEGER :: IERR,IPRT,nHit CALL getLine(Track,X0,Phi) CALL countHit(iRowT,Track,nHit,Xrow) IF( nHit == 0 )THEN reso = -1000.0 RETURN ELSEIF( nHit == 1 )THEN ! fit will fail, set position to center of pad reso = Xrow - X0 RETURN ENDIF PCH(1) = "X0 " PCH(2) = "PHI " PCH(3) = "INVRAD " PVAL(1) = X0 PVAL(2) = tan(Phi) PVAL(3) = 0.0 PERR(1) = 0.1 PERR(2) = -0.03 PERR(3) = -0.03 ! normaly no bounds BND1 = 0.0 BND2 = 0.0 ! define the ZTrack ZTrackFit = ZTrack ! define the rows for the fit ALLOCATE( iRowFit(size(iRowT)),STAT=IERR ) iRowFit = iRowT WRITE(UNIT=CHTITL,FMT=*) "+++++ TPC resolution fit, rows:",iRowFit IPRT = 0 CALL MINUIT(TRKCHIZ,CHTITL,NFITP,PCH,PVAL,PERR,PMAT, & BND1,BND2,CHI2,IPRT,IERR) DEALLOCATE(iRowFit) IF( IERR == 0 .OR. IERR == 1 )THEN reso = PVAL(1) - X0 ELSE WRITE(UNIT=IOErr,FMT=*) "WARNING,TRKFIT > resolution fit failed:" & ,IERR," Event:",EventNumber reso = 1000.0 ENDIF END SUBROUTINE TrackFitXLineZ SUBROUTINE TrackFitXLineW(iRowT,Width,Track,reso) INTEGER,DIMENSION(:),INTENT(IN) :: iRowT REAL(KIND=spPrec),INTENT(IN) :: Width TYPE(Line_type),INTENT(IN) :: Track REAL,INTENT(OUT) :: reso CHARACTER(LEN=50) :: CHTITL INTEGER,PARAMETER :: NFITP=4 CHARACTER(LEN=10),DIMENSION(NFITP) :: PCH REAL(KIND=minuitPrec),DIMENSION(NFITP) :: PVAL,PERR,BND1,BND2 REAL(KIND=minuitPrec),DIMENSION(NFITP,NFITP) :: PMAT REAL(KIND=minuitPrec) :: CHI2 REAL(KIND=spPrec) :: X0,Phi,Xrow INTEGER :: IERR,IPRT,nHit CALL getLine(Track,X0,Phi) CALL countHit(iRowT,Track,nHit,Xrow) IF( nHit == 0 )THEN reso = -1000.0 RETURN ELSEIF( nHit == 1 )THEN ! fit will fail, set position to center of pad reso = Xrow - X0 RETURN ENDIF ! is Trackin good? IF( width > 999.0 )THEN reso = -1000.0 RETURN ENDIF PCH(1) = "X0 " PCH(2) = "PHI " PCH(3) = "INVRAD " PCH(4) = "SIGMA " PVAL(1) = X0 PVAL(2) = tan(Phi) PVAL(3) = 0.0 PVAL(4) = Width PERR(1) = 0.1 PERR(2) = -0.03 PERR(3) = -1.0 PERR(4) = -0.3 ! normaly no bounds BND1 = 0.0 BND2 = 0.0 ! define the ZTrack ZTrackFit = Line(0.0,0.0) ! define the rows for the fit ALLOCATE( iRowFit(size(iRowT)),STAT=IERR ) iRowFit = iRowT WRITE(UNIT=CHTITL,FMT=*) "+++++ TPC resolution fit, rows:",iRowFit IPRT = 0 CALL MINUIT(TRKCHI,CHTITL,NFITP,PCH,PVAL,PERR,PMAT, & BND1,BND2,CHI2,IPRT,IERR) DEALLOCATE(iRowFit) IF( IERR == 0 .OR. IERR == 1 )THEN reso = PVAL(1) - X0 ELSE WRITE(UNIT=IOErr,FMT=*) "WARNING,TRKFIT > resolution fit failed:" & ,IERR," Event:",EventNumber reso = 1000.0 ENDIF END SUBROUTINE TrackFitXLineW SUBROUTINE countHit(iRowT,Track,nHit,Xrow) INTEGER,DIMENSION(:),INTENT(IN) :: iRowT TYPE(Line_type),INTENT(IN) :: Track INTEGER,INTENT(OUT) :: nHit REAL(KIND=spPrec),INTENT(OUT) :: Xrow INTEGER :: iPadNear,iP,iPad,iG,iGroup,iGhit,iR,iRow REAL :: Amplitude,T0 REAL(KIND=spPrec) :: Phi,X0 TYPE(Line_type) :: RowTrack CALL getLine(Track,X0,Phi) Xrow = -1000.0 ! is there a hit in the row(s)? nHit = 0 DO iR=1,size(iRowT) iRow = iRowT(iR) DO iG = 1,nGroupsinRow(iRow) iGroup = iGroupofRow(iG,iRow) CALL GetSignalGroup(iGroup,Amplitude,T0) IF( Amplitude > 0.0 )THEN nHit = nHit + 1 iGhit = iGroup ENDIF ENDDO ENDDO IF( nHit == 1 )THEN ! fit will fail, set position to center of pad iPadNear = -1 DO iP=1,nPadsinGroup(iGhit) iPad = iPadofGroup(iP,iGhit) ! this works only for demultiplexed data CALL GetSignal(iPad,Amplitude,T0) IF( Amplitude > 0.0 )THEN iPadNear = iPad ENDIF ENDDO IF( iPadNear>0 )THEN RowTrack = Line(getCentre(PadLocation(iPadNear)),Phi) CALL getLine(RowTrack,Xrow,Phi) ENDIF ENDIF END SUBROUTINE countHit ! ************************************************************* ! calculate the chisq ***************************************** ! ************************************************************* SUBROUTINE TRKCHI(NPAR,GRAD,FVAL,XVAL,IFLAG) ! ------------------------------------------------------------- ! subroutine to be called by MINUIT - calculates the chisq ! as a function of 4 parameters ! X0: X(Y=0) ! PHI: PHI(Y=0) ! INVRAD: 1/RADIUS ! WIDTH: sigma of track ! ------------------------------------------------------------- ! global variables used: ! iRowFit : which rows to fit ! ------------------------------------------------------------- INTEGER,INTENT(IN) :: IFLAG , NPAR REAL(KIND=minuitPrec),DIMENSION(4),INTENT(IN) :: XVAL REAL(KIND=minuitPrec),DIMENSION(4),INTENT(OUT) :: GRAD REAL(KIND=minuitPrec),INTENT(OUT) :: FVAL REAL(KIND=minuitPrec) :: CHISQ REAL(KIND=spPrec) :: X0,Phi,invRad,width,Yrow TYPE(XYTrack_type) :: fitTrack TYPE(Line_type) :: localTrack INTEGER :: iR,iRow,nFitRows,myPar ! dummy to avoid compiler warnings IF( IFLAG==9 ) myPar=NPAR GRAD = 0.0_minuitPrec CHISQ = 0.0_minuitPrec X0 = XVAL(1) Phi = atan(XVAL(2)) invRad = XVAL(3) width = XVAL(4) fitTrack = XYTrack(Line(X0,Phi),invRad) nFitRows = size(iRowFit) IF( nFitRows == 0 )THEN WRITE(UNIT=IOErr,FMT=*) "ERROR, TRKCHI > no rows for fit" FVAL = 1.0E6 RETURN ENDIF DO iR=1,nFitRows iRow = iRowFit(iR) ! calculate expected probability distribution Yrow = Length(getCentre(PadLocation(iPadofRow(1,iRow)))) localTrack = Line(fitTrack,Yrow) SELECT CASE (FitType) CASE (1) CHISQ = CHISQ + ROWCHI(localTrack,width,iRow) CASE (2) CHISQ = CHISQ + ROWCHItrunc(localTrack,width,iRow) CASE (3) CHISQ = CHISQ + ROWCHIPRF(localTrack,width,iRow) CASE DEFAULT PRINT*," ERROR, TRKCHI > No Fit Type defined for",FitType STOP END SELECT ENDDO FVAL = CHISQ END SUBROUTINE TRKCHI SUBROUTINE TRKCHIZ(NPAR,GRAD,FVAL,XVAL,IFLAG) ! ------------------------------------------------------------- ! subroutine to be called by MINUIT - calculates the chisq ! as a function of 3 parameters ! X0: X(Y=0) ! PHI: PHI(Y=0) ! INVRAD: 1/RADIUS ! the width of the track is calculated from the drift distance. ! ------------------------------------------------------------- ! global variables used: ! ZTrackFit : line for Z(Y) used for track width ! iRowFit : which rows to fit ! ------------------------------------------------------------- INTEGER,INTENT(IN) :: IFLAG , NPAR REAL(KIND=minuitPrec),DIMENSION(3),INTENT(IN) :: XVAL REAL(KIND=minuitPrec),DIMENSION(3),INTENT(OUT) :: GRAD REAL(KIND=minuitPrec),INTENT(OUT) :: FVAL REAL(KIND=minuitPrec) :: CHISQ REAL(KIND=spPrec) :: X0,Phi,invRad,width,Yrow TYPE(XYTrack_type) :: fitTrack TYPE(Line_type) :: localTrack INTEGER :: iR,iRow,nFitRows,myPar REAL :: Z ! dummy to avoid compiler warnings IF( IFLAG==9 ) myPar=NPAR GRAD = 0.0_minuitPrec CHISQ = 0.0_minuitPrec X0 = XVAL(1) Phi = atan(XVAL(2)) invRad = XVAL(3) fitTrack = XYTrack(Line(X0,Phi),invRad) nFitRows = size(iRowFit) IF( nFitRows == 0 )THEN WRITE(UNIT=IOErr,FMT=*) "ERROR,TRKCHIZ > no rows for fit" FVAL = 1.0E6 RETURN ENDIF DO iR=1,nFitRows iRow = iRowFit(iR) ! calculate expected probability distribution Yrow = Length(getCentre(PadLocation(iPadofRow(1,iRow)))) Z = getLineX(ZTrackFit,Yrow) width = trackWidthZ(Z) localTrack = Line(fitTrack,Yrow) SELECT CASE (FitType) CASE (1) CHISQ = CHISQ + ROWCHI(localTrack,width,iRow) CASE (2) CHISQ = CHISQ + ROWCHItrunc(localTrack,width,iRow) CASE (3) CHISQ = CHISQ + ROWCHIPRF(localTrack,width,iRow) CASE DEFAULT PRINT*," ERROR, TRKCHIZ > No Fit Type defined for",FitType STOP END SELECT ENDDO FVAL = CHISQ END SUBROUTINE TRKCHIZ FUNCTION ROWCHI(fitTrack,width,iRow) RESULT(CHISQ) ! ------------------------------------------------------------- ! calculates the chisq in one row ! ------------------------------------------------------------- TYPE(Line_type),INTENT(IN) :: fitTrack REAL(KIND=spPrec),INTENT(IN) :: width INTEGER,INTENT(IN) :: iRow REAL(KIND=minuitPrec) :: CHISQ,totAmpl,GinR INTEGER :: iG,iGroup,iP,iPad REAL(KIND=minuitPrec),DIMENSION(maxGinRow) :: expectAmpl REAL :: Amplitude,T0 CHISQ = 0.0_minuitPrec ! calculate expected probability distribution totAmpl = 0.0 expectAmpl = 0.0 DO iG=1,nGroupsinRow(iRow) iGroup = iGroupofRow(iG,iRow) DO iP=1,nPadsinGroup(iGroup) iPad = iPadofGroup(iP,iGroup) ! only pads in this row IF( PadRow(iPad) /= iRow ) cycle ! don't count if demultiplexed CALL getSignal(iPad,Amplitude,T0) IF( Amplitude < 0.0 ) cycle expectAmpl(iG) = expectAmpl(iG) + PRF(fitTrack,width,iPad) ENDDO totAmpl = totAmpl + expectAmpl(iG) ENDDO IF( totAmpl >= 1.0E-8 ) expectAmpl = expectAmpl / totAmpl GinR = nGroupsinRow(iRow) ! calculate contribution to chisq DO iG=1,nGroupsinRow(iRow) iGroup = iGroupofRow(iG,iRow) CALL GetSignalGroup(iGroup,Amplitude,T0) IF( Amplitude > 0.0 )THEN CHISQ = CHISQ & - REAL(Amplitude,minuitPrec)*REAL(ADC2Electron,minuitPrec) & * log( ( expectAmpl(iG) + REAL(pNoise,minuitPrec) + 1.0E-8_minuitPrec ) & / (1.0_minuitPrec + GinR*REAL(pNoise,minuitPrec) ) ) ENDIF ENDDO END FUNCTION ROWCHI FUNCTION ROWCHItrunc(fitTrack,width,iRow) RESULT(CHISQ) ! ------------------------------------------------------------- ! calculates the chisq in one row ! use only the 'nUsePad' pads closest to track ! ------------------------------------------------------------- TYPE(Line_type),INTENT(IN) :: fitTrack REAL(KIND=spPrec),INTENT(IN) :: width INTEGER,INTENT(IN) :: iRow REAL(KIND=minuitPrec) :: CHISQ,totAmpl INTEGER :: iP,iPad,i,j REAL :: Amplitude,T0,dist INTEGER,PARAMETER :: nUsePad=4 REAL(KIND=minuitPrec),DIMENSION(nUsePad) :: expectAmpl REAL,DIMENSION(nUsePad) :: minDist INTEGER,DIMENSION(nUsePad) :: iUsePad CHISQ = 0.0_minuitPrec ! which pads to use in this row? minDist = 9999.0 iUsePad = 0 DO iP=1,nPadsinRow(iRow) iPad = iPadofRow(iP,iRow) dist = Distance( PadLocation(iPad) , fitTrack ) whichPad : DO i=1,nUsePad IF( dist pad not used" cycle ENDIF ! don't count if demultiplexed CALL getSignal(iPad,Amplitude,T0) IF( Amplitude < 0.0 ) cycle expectAmpl(iP) = PRF(fitTrack,width,iPad) totAmpl = totAmpl + expectAmpl(iP) ENDDO IF( totAmpl >= 1.0E-8 ) expectAmpl = expectAmpl / totAmpl ! calculate contribution to chisq DO iP=1,nUsePad iPad = iUsePad(iP) IF( iPad==0 ) cycle CALL getSignal(iPad,Amplitude,T0) IF( Amplitude < 0.0 ) cycle CHISQ = CHISQ & - REAL(Amplitude,minuitPrec)*REAL(ADC2Electron,minuitPrec) & * log( ( expectAmpl(iP) + REAL(pNoise,minuitPrec) + 1.0E-8_minuitPrec ) & / (1.0_minuitPrec + REAL(nUsePad)*REAL(pNoise,minuitPrec) ) ) ENDDO END FUNCTION ROWCHItrunc FUNCTION ROWCHIPRF(fitTrack,width,iRow) RESULT(CHISQ) ! ------------------------------------------------------------- ! calculates the chisq in one row assuming error(A) ~ Amplitude ! ------------------------------------------------------------- ! WARNING: for multiplexed data ! use only truncated version or deMultiplex ! ------------------------------------------------------------- TYPE(Line_type),INTENT(IN) :: fitTrack REAL(KIND=spPrec),INTENT(IN) :: width INTEGER,INTENT(IN) :: iRow REAL(KIND=minuitPrec) :: CHISQ,CHI INTEGER :: iGmax,iP,iPad,iUse,Middle,Left,Right,nUsePad REAL :: Amplitude,T0,Fac,SumA,SumB,minD,Dist,ALeft,ARight INTEGER,DIMENSION(:),ALLOCATABLE :: iUsePad ! set nUsePad < maxPinRow for truncation ! nUsePad = maxPinRow -> no truncation nUsePad = maxPinRow ! nUsePad = 5 ALLOCATE( iUsePad(nUsePad) ) CHISQ = 0.0_minuitPrec IF( nUsePad == maxPinRow )THEN iUse = nPadsinRow(iRow) iUsePad = iPadofRow(:,iRow) ELSE CALL GmaxSignalRow(iRow,iGmax,Amplitude,T0) iUse = 0 Middle = 0 minD = 1.0E6 DO iP=1,nPadsinGroup(iGmax) iPad = iPadofGroup(iP,iGmax) IF( PadRow(iPad) /= iRow) cycle Dist = Distance(PadLocation(iPad),fitTrack) IF( Dist < minD )THEN minD = Dist Middle = iPad ENDIF ENDDO IF( Middle/=0 ) THEN iUse = iUse+1 iUsePad(iUse) = Middle CALL Neighbour(Middle,Left,Right) CALL getSignal(Left,ALeft,T0) CALL getSignal(Right,ARight,T0) DO IF( iUse>=nUsePad )EXIT IF( ALeft<=0.0 .AND. ARight<=0.0 )EXIT IF( ALeft > ARight )THEN Middle = Left iUse = iUse+1 iUsePad(iUse) = Middle CALL Neighbour(Middle,Left,Right) CALL getSignal(Left,ALeft,T0) ELSE Middle = Right iUse = iUse+1 iUsePad(iUse) = Middle CALL Neighbour(Middle,Left,Right) CALL getSignal(Right,ARight,T0) ENDIF ENDDO ENDIF ENDIF ! normalisation SumA = 0.0 SumB = 0.0 DO iP=1,iUse iPad = iUsePad(iP) CALL getSignal(iPad,Amplitude,T0) ! don't count if demultiplexed IF( Amplitude < 0.0 ) cycle SumA = SumA + Amplitude*PRF(fitTrack,width,iPad)/AERR(Amplitude)**2 SumB = SumB + (Amplitude/AERR(Amplitude))**2 ENDDO IF( iUse<1 .OR. SumB<=0.0 )THEN CHISQ = 0.0 DEALLOCATE( iUsePad ) RETURN ENDIF Fac = SumA/SumB IF( Fac<=1.0E-7 )THEN CHISQ = 1.0E13 DEALLOCATE( iUsePad ) RETURN ENDIF ! calculate contribution to chisq DO iP=1,iUse iPad = iUsePad(iP) CALL GetSignal(iPad,Amplitude,T0) ! don't count if demultiplexed IF( Amplitude < 0.0 ) cycle CHI = ( Amplitude - PRF(fitTrack,width,iPad)/Fac )/AERR(Amplitude) CHISQ = CHISQ + CHI**2 ENDDO DEALLOCATE( iUsePad ) END FUNCTION ROWCHIPRF FUNCTION AERR(A) RESULT(ERROR) ! this is a guestimate of the error on the amplitude REAL,INTENT(IN) :: A REAL :: ERROR IF( A<=0.0 )THEN ERROR = minSignal ELSE IF( size(AERRpar) > 0 )THEN ! take the first value as constant term ERROR = POLY(A,AERRpar(2:size(AERRpar))) ERROR = SQRT( ERROR**2 + AERRpar(1)**2 ) IF( ERROR<=0.0 ) ERROR = A*0.1 ELSE ERROR = A*0.1 ENDIF ENDIF END FUNCTION AERR FUNCTION PRF(track,Width,iPad) RESULT(integral) TYPE(Line_type),INTENT(IN) :: track REAL(KIND=spPrec),INTENT(IN) :: Width INTEGER,INTENT(IN) :: iPad INTEGER :: iR REAL(KIND=minuitPrec) :: integral REAL(KIND=spPrec) :: X0,Phi,bias,XTrack TYPE(Coordinate_type) :: point TYPE(Line_type) :: fittrack ! bias correction iR = iRowBias(PadRow(iPad)) IF( iR /= 0 )THEN CALL getLine(track,X0,Phi) point = getCentre(PadLocation(iPad)) XTrack = getLineX(track,Length(point)) bias = SPLINE(biasX,biasData(:,iR),XTrack) fittrack = Line(X0+bias,Phi) ELSE fittrack = track ENDIF IF( LPRF )THEN integral = QRAT(fitTrack,width,iPad) ELSE integral = intTrackPad(fitTrack,width,PadLocation(iPad)) ENDIF END FUNCTION PRF FUNCTION QRAT(track,width,iPad) RESULT(integral) ! PadResponse function as determined from data ! Width of track is FWHM TYPE(Line_type),INTENT(IN) :: track REAL(KIND=spPrec),INTENT(IN) :: width INTEGER,INTENT(IN) :: iPad TYPE(Coordinate_type) :: point REAL(KIND=minuitPrec) :: integral REAL(KIND=spPrec) :: y,x,c,gam,del,a,b,a2,a4,b2,b4,p,q,qn,qd REAL :: Z point = getCentre(PadLocation(iPad)) ! if Ztrack is provided for the fit this should be reasonable ! otherwise Z should be 0 Z = getLineX(ZTrackFit,Length(point)) del = POLY(Z,PRFdel) a = PRFa b = PRFb c = DISTANCE(track,point) gam = Width IF( del < 0.0 )THEN PRINT*,"ERROR QRAT> parameters not provided" integral = 0.0 RETURN ENDIF IF( c>del*0.5 )THEN integral = 0.0 RETURN ENDIF y = gam / 2.0 x = del / 2.0 a4 = a / x**4 p =-(1.0+a) a2 = p / x**2 b4 = b / y**4 q = 1.0 + 2.0*a2*y**2 + (2.0*a4-b4)*y**4 b2 = q / y**2 qn = 1.0 + a2*c**2 + a4*c**4 qd = 1.0 + b2*c**2 + b4*c**4 integral = qn / qd if(integral<0.0 )then PRINT*, "WARNING: QRAT<0",c,gam endif END FUNCTION QRAT ! <<<< MINUIT <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< END MODULE TPCtrack