MODULE MYanalyzer USE IOunits USE Threshold USE TPCdata USE TPCpaw USE TPCcalcAmplitude USE TPCtrack USE MyMinuit USE Space PRIVATE PUBLIC :: DenseData PUBLIC :: analyzEvent,MY_book PRIVATE :: NF_analysis,HF_test,HF_over PRIVATE :: noOverflow,determinePRF,fitAmplRow,PRFCHI CONTAINS SUBROUTINE analyzEvent(iEvent) INTEGER,INTENT(IN) :: iEvent INTEGER :: IERR,IVETO,nOver REAL,DIMENSION(100) :: Xdata TYPE(Line_type) :: seedTrack,ZTrack,fitTrack REAL(KIND=spPrec) :: Z,THETA INTEGER,DIMENSION(10) :: iGroupV=(/8,14,20,26,32,34,40,46,52,58/) INTEGER,DIMENSION(1) :: iRowSeed1=(/2/),iRowSeed2=(/6/) INTEGER,DIMENSION(5) :: iRowT=(/2,3,4,5,6/) ! Where am I? ! PRINT*,"Event: ",iEvent,EventNumber," Run:",RunNumber ! get the signal amplitudes CALL calcAmplitudeFoil() ! CALL HF_amplitude(iRowT) ! preselection CALL eventVeto(iGroupV,iRowSeed1,iRowSeed2,iRowT,IVETO) CALL HF_veto(IVETO) IF( IVETO /= 0 )RETURN nOver = noOverflow(iRowT) CALL HF_over(nOver) IF( nOver /= 0 )RETURN ! plain event dispay and plot pulse shapes ! CALL NF_event() ! CALL plotPulseRow(iRowT) ! CALL plotPulseEvent() ! call hf_test(iEvent) ! find seed track and NF for event display CALL findSeedTrack(iRowSeed1,iRowSeed2,iRowT,seedTrack,IERR) CALL HF_seedTrack(IERR+1) ! CALL NF_event(seedTrack) IF( IERR /= 0 ) RETURN ! Z-Y track, linear regression or fit CALL TrackFitZ(iRowT,ZTrack) CALL getLine(ZTrack,Z,THETA) fitTrack = seedTrack ! X-Y track, with or without variable track width, NF for event display ! fit 3 parameters ! CALL TrackFitWidth(iRowT,fitTrack,Width) ! fit 2 parameters, calculate width from time0 CALL TrackFitXY(iRowT,ZTrack,fitTrack) ! CALL NF_event(fitTrack,Width) IF( .NOT.(fitTrack == seedTrack) ) & CALL determinePRF(fitTrack,Z) END SUBROUTINE analyzEvent FUNCTION noOverflow(iRowT) RESULT(nOver) INTEGER,DIMENSION(:),INTENT(IN) :: iRowT INTEGER :: nOver,iR,iRow,iG,iGroup nOver = 0 ! reject events with overflow DO iR=1,size(iRowT) iRow = iRowT(iR) DO iG=1,nGroupsinRow(iRow) iGroup = iGroupofRow(iG,iRow) IF( minval(IADC(1:nTbGroup(iGroup),iGroup)) <= 0 ) nOver = nOver+1 ENDDO ENDDO END FUNCTION noOverflow SUBROUTINE determinePRF(fitTrack,Z) TYPE(Line_type),INTENT(IN) :: fitTrack REAL(KIND=spPrec),INTENT(IN) :: Z INTEGER :: iR,iRow,iP,iPad,iGmax REAL :: Amplitude,amplrow,XPad,Ypad,xTrack,deltaX,deltaT REAL :: T0,Trow,AmplMax REAL,DIMENSION(100) :: Xdata INTEGER,DIMENSION(3) :: iRowTest=(/3,4,5/) TYPE(Coordinate_type) :: Point REAL(KIND=spPrec) :: X0Fit,PhiFit CALL getLine(fitTrack,X0Fit,PhiFit) DO iR = 1,size(iRowTest) iRow = iRowTest(iR) CALL GmaxSignalRow(iRow,iGmax,AmplMax,Trow) IF( AmplMax>200.0 ) CYCLE CALL GetSignalRow(iRow,AmplRow,T0) IF( AmplRow<30.0 ) CYCLE CALL fitAmplRow(fitTrack,Z,iRow,AmplRow) DO iP = 1,nPadsinRow(iRow) iPad = iPadofRow(iP,iRow) CALL GetSignal(iPad,Amplitude,T0) point = getCentre(PadLocation(iPad)) Xpad = Pitch(Point) Ypad = Length(Point) Xtrack = getLineX(fitTrack,Ypad) deltaX = Distance(fitTrack,Point) IF( Xpad 0.001 .OR. ABS(deltaX)<6.0 )THEN XDATA(1) = X0Fit XDATA(2) = PhiFit XDATA(3) = Z XDATA(4) = deltaX XDATA(5) = deltaT XDATA(6) = Amplitude/AmplRow XDATA(7) = AmplRow XDATA(8) = iPad CALL HFN(952,Xdata) ENDIF ENDDO ENDDO END SUBROUTINE determinePRF SUBROUTINE fitAmplRow(fitTrack,Z,iRow,AmplRow) TYPE(Line_type),INTENT(IN) :: fitTrack REAL(KIND=spPrec),INTENT(IN) :: Z INTEGER,INTENT(IN) :: iRow REAL,INTENT(OUT) :: AmplRow INTEGER :: IPRT,IERR REAL(KIND=spPrec) :: X0,Phi CHARACTER(LEN=50) :: CHTITL INTEGER,PARAMETER :: NFITP=5 CHARACTER(LEN=10),DIMENSION(NFITP) :: PCH REAL(KIND=minuitPrec) :: CHISQ REAL(KIND=minuitPrec),DIMENSION(NFITP) :: PVAL,PERR,BND1,BND2 REAL(KIND=minuitPrec),DIMENSION(NFITP,NFITP) :: PMAT CALL getLine(fitTrack,X0,Phi) PCH(1) = "AmplRow " PCH(2) = "X0 " PCH(3) = "Phi " PCH(4) = "Sigma " PCH(5) = "Row " PVAL(1) = 30.0 PVAL(2) = X0 PVAL(3) = Phi PVAL(4) = trackWidth(Line(Z,0.0_spPrec),0.0) PVAL(5) = iRow PERR(1) = 3.0 PERR(2) = -0.1 PERR(3) = -1.0 PERR(4) = -1.0 PERR(5) = -1.0 ! normaly no bounds BND1 = 0.0 BND2 = 0.0 WRITE(UNIT=CHTITL,FMT=*) "+++++ PRF fit, row:",iRow IPRT = 0 CALL MINUIT(PRFCHI,CHTITL,NFITP,PCH,PVAL,PERR,PMAT, & BND1,BND2,CHISQ,IPRT,IERR) AmplRow = PVAL(1) END SUBROUTINE fitAmplRow SUBROUTINE PRFCHI(NPAR,GRAD,FVAL,XVAL,IFLAG) INTEGER,INTENT(IN) :: IFLAG , NPAR REAL(KIND=minuitPrec),DIMENSION(5),INTENT(IN) :: XVAL REAL(KIND=minuitPrec),DIMENSION(5),INTENT(OUT) :: GRAD REAL(KIND=minuitPrec),INTENT(OUT) :: FVAL REAL(KIND=minuitPrec) :: CHISQ,sum TYPE(Line_type) :: Track INTEGER :: iRow,iP,iPad,myPar REAL(KIND=spPrec) :: X0,Phi,width REAL :: AmplRow,Amplitude,T0 ! dummy to avoid compiler warnings IF( IFLAG==9 ) myPar=NPAR GRAD = 0.0_minuitPrec CHISQ = 0.0_minuitPrec AmplRow = XVAL(1) X0 = XVAL(2) Phi = XVAL(3) width = XVAL(4) iRow = XVAL(5) Track = Line(X0,Phi) DO iP = 1,nPadsinRow(iRow) iPad = iPadofRow(iP,iRow) CALL GetSignal(iPad,Amplitude,T0) sum = (Amplitude-PRF(Track,width,iPad)*AmplRow)/Aerr(Amplitude) ! sum = 1.0_minuitPrec - PRF(Track,width,iPad)*AmplRow/Amplitude CHISQ = CHISQ + sum**2 ENDDO FVAL = CHISQ END SUBROUTINE PRFCHI SUBROUTINE MY_book() ! book histograms & Ntuples for analyzEvent CHARACTER(LEN=8),DIMENSION(1000) :: CHtag CHARACTER(LEN=20) :: CHtitle INTEGER :: ntag,i ! Histograms, numbers >1039 CALL HBOOK1(1040,"test",64,0.5,64.5,0.0) call hbook1(1041,"t0",400,0.0,4000.0,0.0) call hbook1(1050,"nOver",40,-0.5,39.5,0.0) ! do i=1,ngroups ! WRITE(UNIT=CHtitle,FMT=*) "mean amplitude",i ! CALL HBPROF(1100+i,CHtitle,50,0.0,10000.0,-0.5,4000.0,' ') ! enddo CALL HBOOK2(1046,"time difference",40,0.0,0.5,40,-50.0,350.0,0.0) CALL HBOOK2(1047,"shape",40,0.0,0.5,40,0.0,1.0,0.0) CALL HBOOK2(1048,"shape",40,-50.0,350.0,40,0.0,1.0,0.0) CHtag(1) = "X0FIT" CHtag(2) = "PHIFIT" CHtag(3) = "Z" CHtag(4) = "DELTAX" CHtag(5) = "DELTAT" CHtag(6) = "AMPLNORM" CHtag(7) = "AMPLROW" CHtag(8) = "PAD" ntag = 8 CALL NT_book(952,"PRF",ntag,CHtag) END SUBROUTINE MY_book ! NTUPLE fill ********************************************** SUBROUTINE NF_analysis(Xdata) REAL,DIMENSION(:),INTENT(IN) :: Xdata CALL HFN(951,Xdata) END SUBROUTINE NF_analysis ! HISTOGRAM fill ********************************************** SUBROUTINE HF_test(iEvent) INTEGER,INTENT(IN) :: iEvent integer:: iRow,iG,iGroup,nHitRow,iTmax,iT0 real :: amplitude,t0,amplsum,amplratio,ampllate,amplmax DO iGroup=1,nGroups CALL getSignalGroup(iGroup,Amplitude,T0) if( amplitude>1.0 )then CALL HFILL(1040,real(iGroup),0.0,amplitude) call hfill(1041,t0,0.0,1.0) endif ! call hfill(1100+iGroup,REAL(iEvent),amplitude,1.0) enddo DO iRow=1,nRows amplsum = 0.0 amplmax = 0.0 nHitRow = 0 DO iG=1,nGroupsinRow(iRow) iGroup = iGroupofRow(iG,iRow) CALL getSignalGroup(iGroup,Amplitude,T0) amplsum = amplsum + Amplitude IF( Amplitude > amplmax )THEN iTmax = T0 amplmax = Amplitude ENDIF IF( Amplitude > 0.0 ) nHitRow = nHitRow + 1 ENDDO IF( nHitRow > 1 )THEN DO iG=1,nGroupsinRow(iRow) iGroup = iGroupofRow(iG,iRow) CALL getSignalGroup(iGroup,Amplitude,T0) IF( Amplitude > 0.0 )THEN amplratio = Amplitude / amplsum iT0 = T0 IF( nTbin/=0 )THEN ampllate = sum(ADC(iT0+81:iT0+130,iGroup)) / 50.0 ampllate = -ampllate / Amplitude ELSE ampllate = 0.0 ENDIF CALL HF2(1046,amplratio,REAL(iT0-iTmax),1.0) CALL HF2(1047,amplratio,ampllate,1.0) CALL HF2(1048,REAL(iT0-iTmax),ampllate,1.0) ENDIF ENDDO ENDIF ENDDO END SUBROUTINE HF_test SUBROUTINE HF_over(nOver) INTEGER,INTENT(IN) :: nOver CALL HF1(1050,REAL(nOver),1.0) END SUBROUTINE HF_over ! >>>> WRITE DENSE DATA >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> SUBROUTINE DenseData() USE DataFormats, only : WriteDenseData ! INTEGER :: IVETO ! get the signal amplitudes CALL calcAmplitudeFoil() ! preselection ! CALL eventVeto(IVETO) ! IF( IVETO /= 0 )RETURN CALL WriteDenseData() END SUBROUTINE DenseData END MODULE MYanalyzer