*********************************************** * * SUBROUTINE MDCDIGI * * * * * last updated on: 28/11/2008 by I. Koenig * *********************************************** IMPLICIT NONE #include "drtups.inc" #include "geant321/gcflag.inc" #include "user.inc" #include "kinetups.inc" INTEGER*4 NVDIM, NHDIM, NHMAX, NSEC, NMOD, NLAY, NTR, NH, NTL PARAMETER (NVDIM=4, NHDIM=16, NHMAX=600, NSEC=6, NMOD=4, NLAY=6+1) CHARACTER*4 IUDET INTEGER*4 ITRACK(NHMAX), NUMBV(NVDIM,NHMAX), NHITS REAL*4 HITS(NHDIM,NHMAX) INTEGER*4 NUMVS(NVDIM) INTEGER*4 JSEC, JMOD, JLAY, JHIT, JTRK, I, J, K, L INTEGER*4 TRK(MAXHIT) INTEGER*4 TRKHIT(MAXHIT) c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ INTEGER*4 IN_SECTOR, N_PLANES_HIT(4) REAL*4 PTOT, THETA LOGICAL IS_HIT, IS_ACCEPTED c--------------------------------------------------------------------------- * extract drift-chamber hits from 'JHITS' c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ IF(IMACC.EQ.1) THEN ! do some preparation acceptance-matrix calculation DO I=1,4 N_PLANES_HIT(I) = 0 ENDDO IN_SECTOR = INT(PHI_EVENT/60.0) ! find out which sector is hit IF(IN_SECTOR.EQ.0) IN_SECTOR = 6 ENDIF c--------------------------------------------------------------------------- DO 20 JSEC=1,NSEC ! loop over sectors (=6) NUMVS(1) = JSEC NUMVS(2) = 0 NUMVS(3) = 0 NUMVS(4) = 0 ! ==> select all hits on sector JSEC in gfhits(...) DO 25 JMOD=1,NMOD ! loop over drift chambers (=4) NTR = 0 NH = 0 DO 30 JLAY=1,NLAY ! loop over wire planes (6x sense + 1 cathode) NHITS = 0 NTL = 0 ! number of hits in a single layer IUDET = 'D'//CHAR(JMOD+48)//'S'//CHAR(JLAY+48) ! ==> select det IF(JLAY.EQ.NLAY) IUDET = 'D'//CHAR(JMOD+48)//'C4' ! 4th cathode CALL GFHITS('DRIF',IUDET,NVDIM,NHDIM,NHMAX,0,NUMVS, & ITRACK,NUMBV,HITS,NHITS) ! retrieve hits from Geant store IF (NHITS.GT.NHMAX) THEN ! overflow WRITE(6,*) 'ERROR: overflow in mdcdigi.F (nHits>', & NHMAX,') in MDC layer',JSEC,JMOD,JLAY RETURN ENDIF IS_HIT = .FALSE. DO 40 JHIT=1,NHITS ! loop over hits in Geant store JTRK = ITRACK(JHIT) IF(JTRK.LE.0) THEN WRITE(6,*) 'ERROR: ITRA=0 in mdcdigi.F' ! this should never happen... RETURN ENDIF DO K=1,NTR ! loop over tracks seen so far (in all layers) IF(TRK(K).EQ.JTRK) THEN TRKHIT(K) = TRKHIT(K)+1 ! increment track hit counter IF(TRKHIT(K).GT.MAXMUL) THEN IF(TRKHIT(K).EQ.MAXMUL+1) THEN ! print message once only WRITE(6,*) 'Track with too many hits (>', & MAXMUL, ') in MDC',JSEC,JMOD ENDIF GOTO 40 ENDIF GOTO 401 ENDIF ENDDO IF(NTR.LT.MAXHIT) THEN NTR = NTR + 1 TRK(NTR) = JTRK ! new track, add to list TRKHIT(NTR) = 1 ENDIF 401 NH = NH + 1 NTL = NTL + 1 IF(NTL.GT.MAXTRKMDC) THEN ! cannot handle this many... WRITE(6,*) 'Too many Tracks in MDC layer', & JSEC,JMOD,JLAY RETURN ENDIF c Here come the goodies that go into the tree c MDCTRAK(NH) = JTRK ! track number LAYER(NH) = JLAY ! wire plane number c MDCX(NH) = LUNIT*HITS(12,JHIT) ! distance along wire c MDCY(NH) = LUNIT*HITS(13,JHIT) ! distance to wire MDCX(NH) = LUNIT*HITS(7,JHIT) ! y distance (in mm) MDCY(NH) = LUNIT*HITS(8,JHIT) ! x distance (in mm) MDCTHETA(NH) = HITS(15,JHIT) ! angle of incidence MDCPHI(NH) = HITS(16,JHIT) ! azimuthal angle MDCTOF(NH) = HITS(11,JHIT)*1.e9 ! tof (in ns) MDCMOM(NH) = 1.e3*SQRT(HITS(4,JHIT)**2+ & HITS(5,JHIT)**2+HITS(6,JHIT)**2) ! momentum (in MeV/c) c DO K=1,NHDIM ! clean-up Geant store c HITS(K,L) = 0.0 c ENDDO c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ IF(JTRK.EQ.1) THEN ! this is the primary track IS_HIT = .TRUE. ENDIF c--------------------------------------------------------------------------- c if (eventid.eq.154) write(6,4111) MDCTRAK(NH),JSEC,JMOD,JLAY, c & MDCTOF(NH), MDCMOM(NH) c 4111 format('Trk:',i2,' (s,m,l)=',3i2,' tof= ',f10.6,' p= ',f10.6) c--------------------------------------------------------------------------- c write(6,*) NH,MDCTRAK(NH),JSEC,JMOD,LAYER(NH),NHITS 40 CONTINUE ! end loop over hits c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ IF(IMACC.EQ.1 .AND. JSEC.EQ.IN_SECTOR .AND. IS_HIT + .AND. JLAY.LT.NLAY) THEN N_PLANES_HIT(JMOD) = N_PLANES_HIT(JMOD)+1 ! count planes hit ENDIF c--------------------------------------------------------------------------- 30 CONTINUE ! end loop over wire planes NHIT = NH IF(IMDC.GE.1) THEN IF(IROOT.EQ.0) THEN CALL HFNT(3000+JSEC*10+JMOD) ! write data to ntuple ELSE CALL FILLMDC(JSEC,JMOD) ENDIF ENDIF 25 CONTINUE ! end loop over drift chambers c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ IF(IMACC.EQ.1 .AND. JSEC.EQ.IN_SECTOR) THEN ! inc. acceptance matrix IS_ACCEPTED = .TRUE. DO JMOD=1,4 ! look if all 4 MDCs are hit C IF(N_PLANES_HIT(JMOD).LT.5) IS_ACCEPTED = .FALSE. IF(N_PLANES_HIT(JMOD).LT.6) IS_ACCEPTED = .FALSE. END DO DO I=1,NTR IF(TRACK(I).EQ.1) GOTO 251 ! locate 1st i.e. primary track END DO 251 PTOT = SQRT(PMOM(1,I)**2+PMOM(2,I)**2+PMOM(3,I)**2) IF(PTOT.GT.0.0) THETA = 57.295*ACOS(PMOM(3,I)/PTOT) IF(IS_ACCEPTED) THEN ! this is a good one CALL HF2(21,PTOT,THETA,1.0) ! count detected particle ENDIF ENDIF c--------------------------------------------------------------------------- 20 CONTINUE ! end loop over sectors RETURN END