c $Id: gbr2c.f,v 1.2 2001/08/01 20:45:07 sanders Exp $ c $Log: gbr2c.f,v $ c Revision 1.2 2001/08/01 20:45:07 sanders c Fixed dimension error of spot array (was 6, change to 9). c c Revision 1.1.1.1 2001/06/12 11:11:42 cholm c Initial import of BRAG (fomerly known as GBRAHMS c c Revision 1.29 2001/05/09 19:14:02 videbaek c Added another loop to be able to include also primary particle in the list. c Can be enabled using saveoption(2) following the init_detector etc. c c Revision 1.28 2000/12/05 17:22:20 videbaek c Modified TOFW to have first panel with 20 and subsequent with 21. A large re-write of several c places in code. c c Revision 1.27 2000/12/02 19:11:15 videbaek c modify to match geant common block size c c Revision 1.26 2000/10/30 15:24:40 videbaek c added run_mrs c c Revision 1.25 2000/05/25 18:04:52 videbaek c Fix error in volume setting for TOFW (iei TFP1). c c Slatnumber too has to be checked. c --------------------------------------------------------------------- c c Revision 1.24 2000/04/03 18:50:45 videbaek c Changed default names c c c Revision 1.23 1999/12/17 19:17:39 videbaek c Insert peoper number of TOFw slats per module c c Revision 1.22 1999/02/01 21:09:56 videbaek c Added mult in hits c c Revision 1.21 1999/02/01 20:57:01 videbaek c added hits for multiplicity c c Revision 1.20 1998/12/01 21:50:21 videbaek c Updated hit information for panel TOFW c c c Revision 1.17 1998/11/12 22:12:21 videbaek c fixed end of event output c c Revision 1.16 1998/11/12 19:58:12 videbaek c Generate and read modified .cdat files format c c Revision 1.15 1998/07/16 21:46:08 videbaek c fix type and statement palcement c c Revision 1.14 1998/07/08 17:38:56 videbaek c Modified common block area ws area for hit storage c c Revision 1.12 1998/03/13 20:46:48 videbaek c Updated for cmodeopen() c c Revision 1.11 1998/01/08 18:47:01 videbaek c Modified volume structure for c-files c c Revision 1.10 1997/12/29 19:59:58 videbaek c updated geometry c c Revision 1.8 1997/11/05 22:01:26 videbaek c matching version (t1pr,t2pr c c SUBROUTINE GBR2C integer fileid common /c_out/ fileid write(6,*) '**************************************' write(6,*) 'GBR2C ENTRIES FOR GBRAHMS' write(6,*) ' Version May 9, 2001' write(6,*) 'event: GBRANA' write(6,*) 'event: init_detector' write(6,*) 'Finish: GBREND' write(6,*) 'Book: GBRFILE(file)' write(6,*) '**************************************' fileid=-1 return end SUBROUTINE GBRANA * * This is the next development in the series of program/routines * which can * save output from the GBRAHMS geant. This protottype is likely close * * * to a version which can save hits in a general way, and which may be * included into the standard GEANT as an output option for hits etc. * It write the hits, track information on a stream i/o C-file. It is not * a platform independent file. * I (FV) have modified SONATA such it can read these files. There is also a * C++ class which can read and save the structures. * The latter class can be found in the sub-directory gbr * * * This analysis program is intended when developed further * to be included into the * GBRAHMS program such that we have an uniform output description. * This output is saved on one file with a general volume, and setup description. * This description can be used either directly or to check a given setup. * * Updates * ======= * December * * November 1998 * Do not output empty records to cut down on file sizes. * Do also add general event information like * vertex, impact parameter to have at least some more overall * diagnostics information. * * October 3 1997: * Add a simple first word (i.e. 1) which can be used * to determine if byte-swapping is needed or not. * July 7, 1998 * Add the geant PID to the Hit structure. This will make * the analysis much easier later, and is needed for detector * simulations. * * May 25, 2000: * Fv. The TOFW volume volume setting was fouled up. The proper logic * was there, but somehow another statement that belong above * ws below. Moved statement. Code now checks out for TPM1, TPM2 TOFW * (i.e. TFP1) * May 9m 2001 * FV. Upgrade track selction to be more along the lines of original toughts * Needed for accpetance calculation to be able to keep all generated as wel * as tracks with hit. * PARAMETER (KWBANK=69000,KWWORK=5200) COMMON/GCBANK/NZEBRA,GVERSN,ZVERSN,IXSTOR,IXDIV,IXCONS,FENDQ(16) + ,LMAIN,LR1,WS(KWBANK) DIMENSION IQ(2),Q(2),LQ(8000),IWS(2) EQUIVALENCE (Q(1),IQ(1),LQ(9)),(LQ(1),LMAIN),(IWS(1),WS(1)) * PARAMETER (NVDIM=20,NHMAX=1500,NHDIM=10) DIMENSION NUMVS(NVDIM),NUMBV(NVDIM,NHMAX),HITS(NHDIM,NHMAX), + ITRA(NHMAX) EQUIVALENCE (WS(5000),NUMVS(1)),(WS(5021),NUMBV(1,1)), + (WS(35021),HITS(1,1)),(WS(50021),ITRA(1)) * * COMMON/GCNUM/NMATE ,NVOLUM,NROTM,NTMED,NTMULT,NTRACK,NPART + ,NSTMAX,NVERTX,NHEAD,NBIT * COMMON/GCFLAG/IDEBUG,IDEMIN,IDEMAX,ITEST,IDRUN,IDEVT,IEORUN + ,IEOTRI,IEVENT,ISWIT(10),IFINIT(20),NEVENT,NRNDM(2) * * * Event generator information * integer maxptl integer nptls, ntry integer nnproj,npproj,naproj integer nntarg,nptarg,natarg real bimevt c c common /zbevent/ nptls, bimevt, & nntarg, nptarg, natarg, & nnproj, npproj, naproj, & ntry * real vertpos(3) common /gbevent/ vertpos * * * Geant common blocks * INTEGER UFLAGS(10),NKINE(10) INTEGER NEVSPC,NEVSPP,NEVUSR,NEVPAR REAL SPOT(9), ROTANG COMMON/GCUSER/UFLAGS,NKINE,SPOT & , NEVSPC,NEVSPP,NEVUSR,NEVPAR,ROTANG * CHARACTER*80 CHKINE,CHSTOR,CHREAD COMMON /GCUSRC/ CHKINE,CHSTOR,CHREAD * INTEGER IPART(2),NVERT,NWBUF REAL PVERT(4),UBUF(10),P1(1:4),P2(1:4) REAL VERT(3),TOFGN,p,psq INTEGER ORGTRK,NGEN INTEGER INDX,ISTAT INTEGER IEVT data ievt/0/ integer ibuf(10) equivalence(ubuf, ibuf) c c integer max_volume parameter (max_volume=50) character*4 detector_names(max_volume,10) character*4 detector_set(max_volume) integer detector_list(max_volume,20) integer detector_nsub(max_volume) integer number_detectors c common /setup/ $ detector_list, detector_names, $ detector_set, $ detector_nsub, number_detectors c logical save_all, save_charged, save_primary common /out_option/ save_all, save_charged, save_primary integer fileid common /c_out/ fileid character *4 name1 CHARACTER*4 NAME2, sublist(20) character *4 setname character *4 detname Character*20 napart * integer debug vector debug(5) vector iplot(1) integer num_tracks integer tracklist common /track_list/ num_tracks, tracklist(500) * integer nlev integer lnam(10), lnum(10), lnum_volu(10) logical first_hit * * structures for saving track and hit information * integer Itrack_s(14) real Rtrack_s(14) equivalence(Itrack_s, Rtrack_s) * integer type_offset integer type_subevent integer type_runnumber integer type_vertexpos integer type_impactpar parameter (type_offset = 1) parameter (eventno_offset = 2) parameter (subevent_offset = 3) parameter (runnumber_offset= 4) parameter (vertexpos_offset= 5) parameter (impactpar_offset= 6) integer Ievt_s(10) real Revt_s(10) equivalence(Ievt_s, Revt_s) * integer Ihits_s(22) real Rhits_s(22) equivalence(Ihits_s, Rhits_s) * * integer cwrite * integer lun,lun_d * * HISTOGRAM THE RESULTS FOR GENERATED TRACKS... * lun = 6 lun_d = 6 if(debug(1) .ge. 1) then print *,'entering ntana', nptls print *, bimevt endif ievt = ievt + 1 c c for now always zero c num_tracks = 0 if(debug(1) .eq. 1) then if(ntrack .ge. 1) then WRITE(6,10001) IEVT, NTMULT, NTRACK,NVERTX 10001 FORMAT(1X,'ANALYSIS : IEVT =',I6,' NTMULT=',I4, & ' NTRACK=',I5,' NVERTX=',I5) end if end if c c Setup the event information c August 4 1998 . The event does not contain the Run Number c This was not too smart. Should have been there- modify after MDCI c first_hit = .true. Ievt_s(type_offset) = 1 Ievt_s(eventno_offset) = Nevusr Ievt_s(subevent_offset) = Nevpar Ievt_s(runnumber_offset) = IDRUN Revt_s(vertexpos_offset) = vertpos(3) Revt_s(impactpar_offset) = bimevt if(abs(nevpar) .eq. 1) then nb = 6*4 nbw = cwrite(fileid, Ievt_s, nb) if(nbw .ne. nb) then write(6,*) 'output error', nb, nbw endif first_hit = .false. endif c c If this is the first sub event of a new event the c event header should be output regard less of c track/hit status. c c c Loop over chanber information c ============================= if(debug(1)) then write(6,*) 'Number detectors', number_detectors end if do ichm = 1, number_detectors isub = detector_nsub(ichm) if(debug(1).ge.2) then write(6,*) 'Detector Number#', ichm write(6,*) 'Number of subdet', isub WRITE(6,10001) IEVT, NTMULT, NTRACK,NVERTX end if do is = 1, isub sublist(is) = detector_names(ichm, is) end do * * GET THE HIT INFORMATION * from chamber ichm * if(ntrack .gt. 0) then do is = 1, isub setname = detector_set(ichm) detname = sublist(is) call zfpath(setname, detname, numvs, nlev, lnam, lnum) do k = 1, nlev lnum_volu(k) = lnum(k) end do if(debug(4).ge.2) then write(lun_d,*) '*** Volume list ***' do k = 1, nlev call uhtoc(lnam(k),4,name1,4) write(lun_d,*) name1, lnum(k) end do endif c c in the case of a non existing detector and/or set (nlev <=0) c the remaining part of this loop is skipped c if(nlev.le.0) goto 997 c c Comment for this choice. c lnum(nlev) is 0 for sub divided detectors. c This is the algorithm to calculate local coordinate not in the c subdivided volume but in the 'master'. c This is certianly the right thing to do in case of a virtually segmented c gas detector (like a TPC). But may not be the right thing for a c TOF wall ?! c if(lnum(nlev).eq.0) then if(debug(1).ge.1) then write(6,*) '*******Nlev', nlev write(6,*) setname,' ', detname, nlev endif if(setname .ne. 'H1 ' & .and. setname .ne. 'H2 ' $ .and. setname .ne. 'TFP2' ) then nlev=nlev-1 else c write(6,*) '*******Nlev', nlev c write(6,*) setname,' ', detname, nlev endif endif ITRS = 0 if(debug(3) .eq. 1) write(6,*) 'Entering gfhitc' call uctoh(detname,iudet,4,4) call uctoh(setname,iuset,4,4) CALL GFHITC (iuset, iudet, ITRS, NHITS, ISTAT) if(istat .ne. 1) then write(6,*) '*** GFHITs array exceeded limits', nhits write(6,*) '*** Debug information *****' write(6,*) setname, detname, nlev WRITE(6,10003) NHITS,ichm, detname, istat 10003 FORMAT(1X,'NUMBER OF HITS=',I3,' in chm',i1, & ' (',a4,' stat=',i3) Nhits = 800 if(idebug) then do igfhits = 1, nhits itrack_num = itra(igfhits) WRITE( 6 ,35) ITrack_num, numbv(1,igfhits), & ( HITS( i, igfhits ), i = 1, 9) enddo endif endif DO IGFHITS = 1, NHITS itrack_num = itra(igfhits) CALL GFKINE(itrack_num, $ VERT,PVERT,IPART(1),NVERT,UBUF,NWBUF) CALL GFPART(IPART(1), NAPART, ITRTYP, AMASS, CHARGE, & TLIFE, UBUF, NWBUF) if(save_all .or. $ (save_charged .and. (charge .ne. 0))) then c c is this track in list? c if not add to list of tracks to deal with c if(.not. IsInTrackList(itrack_num)) then num_tracks = num_tracks+1 tracklist(num_tracks) = itrack_num endif c c look at the hits for this module c IF(DEBUG(2).EQ. 1) then WRITE( 6 ,35) ITrack_num,numbv(1,igfhits), & ( HITS( i, igfhits ), i = 1, 9) 35 FORMAT(1h ,' Track number = ',I3,i3,' HITS : ', & 4G12.5,/,1h ,8(4G12.5,/,1h ) ) write(6,*) 'Test', numbv(2,igfhits) endif ihits_s(1) = 2 c c insert track number and volume number c Ihits_s(2) = itrack_num+(abs(nevpar)-1)*10000 ihits_s(3) = detector_list(ichm, is) if(ichm .eq. 30) then ihits_s(3) = detector_list(21, is) endif c c setup volume information so we can convert to local coordinates c call reset_nlevel c if(isub.eq.1) then Ihits_s(19) = numbv(1, igfhits) else Ihits_s(19) = is endif c c This special construct is for TOFW only which has multiple c panels and slats. Note the 21 is hardwired. c Not gotten from any reasonable parameters in the setup. c This requires additional common blocks available to gbr2c.f c if(lnum(nlev) .eq. 0) then lnum_volu(nlev) = numbv(1, igfhits) endif if(lnum(nlev-1) .eq. 0) then lnum_volu(nlev-1) = numbv(1, igfhits) if(lnum(nlev) .eq. 0) then lnum_volu(nlev) = numbv(2, igfhits) endif endif if(setname(1:3) .eq. 'STR') then lnum_volu(nlev-2) = numbv(1, igfhits) endif c-- with only one short panel should be left out c if(setname .eq. 'TFP1') then c Ihits_s(19) = numbv(2, igfhits) c endif if(setname .eq. 'TFP2') then Ihits_s(19) = numbv(2, igfhits) $ + (numbv(1, igfhits)-1)*21 + 20 endif if(debug(5).eq.2) then write(lun_d,*) '*** Modified Volume list ***' do k = 1, nlev call uhtoc(lnam(k),4,name1,4) write(lun_d,*) name1, lnum(k),lnum_volu(k) end do endif call glvolu(nlev, Lnam, lnum_volu, ier) c call gmtod(hits(1,igfhits), Rhits_s(4), 1) call gmtod(hits(4,igfhits), Rhits_s(7), 1) Do i=1,6 Rhits_s(9+i) = hits(i,igfhits) end do Rhits_s(16) = HITS(7,IGFHITS) Rhits_s(17) = HITS(8,IGFHITS) psq = hits(9, igfhits)**2 - amass**2 if(psq.gt. 0) then psq = sqrt(psq) else psq = 0.0 endif Rhits_s(18) = psq Ihits_s(20) = ipart(1) if(first_hit) then nb = 6*4 nbw = cwrite(fileid, Ievt_s, nb) if(nbw .ne. nb) then write(6,*) 'output error', nb, nbw endif first_hit = .false. endif nb = 20*4 nbw = cwrite(fileid, Ihits_s, nb) if(nbw .ne. nb) then write(6,*) 'output error', nb, nbw endif if(debug(4) .eq. 1) then call print_hits(lun, $ setname, detname, ihits_s, rhits_s) end if end if end do 997 end do endif * end loop over chambers * enddo c c Check if additional tracks that are parents of hits C should be included in the primary list c do it = 1, num_tracks indx = tracklist(it) CALL GFKINE(INDX,VERT,PVERT,IPART(1),NVERT,UBUF,NWBUF) CALL GFPART(IPART(1), NAPART, ITRTYP, AMASS, CHARGE, & TLIFE, UBUF, NWBUF) call gfvert(nvert, vert, ntbeam, ngen, tofgn, ubuf, nwbuf) c c There could be two options. c 1) add all decay tracks. c 2) do only add decays tracks if origin was from vertex i.e. c single decay going into a detector hit. c iterate back over tracks.. c if(debug(1) .gt. 0) then write(6,*) 'Iterate..', INDX, ntbeam, ngen endif do while(ntbeam .gt. 0) call uhtoc(ibuf(3),4, name1, 4) call uhtoc(ibuf(4),4, name2, 4) if(name2 .eq. 'DCAY') then indx = ntbeam if(.not. IsInTrackList(ntbeam)) then num_tracks = num_tracks+1 tracklist(num_tracks) = ntbeam if(debug(1) .gt. 0) then write(6,*) 'Add to list', ntbeam endif endif CALL GFKINE(INDX,VERT,PVERT,IPART(1),NVERT,UBUF,NWBUF) CALL GFPART(IPART(1), NAPART, ITRTYP, AMASS, CHARGE, & TLIFE, UBUF, NWBUF) call gfvert(nvert, vert, ntbeam, ngen, tofgn, ubuf, nwbuf) else ntbeam = 0 end if end do end do c c Check if additional tracks that are primaray tracks C but not already in list should be included in the primary list c if(save_primary) then do it = 1, ntrack CALL GFKINE(it,VERT,PVERT,IPART(1),NVERT,UBUF,NWBUF) CALL GFPART(IPART(1), NAPART, ITRTYP, AMASS, CHARGE, & TLIFE, UBUF, NWBUF) call gfvert(nvert, vert, ntbeam, ngen, tofgn, ubuf, nwbuf) if(ntbeam .eq. 0) then if(.not. IsInTrackList(it)) then num_tracks = num_tracks+1 tracklist(num_tracks) = itrack_num if(debug(1) .gt. 0) then write(6,*) 'Add to list #', it endif endif endif end do endif c c c now loop over the tracks which should be included c if(debug(1) .eq. 1) then write(6,*) 'Number of tracks', num_tracks endif c do it = 1, num_tracks c c get the particle and vertex information for this track c ====================================================== indx = tracklist(it) CALL GFKINE(INDX,VERT,PVERT,IPART(1),NVERT,UBUF,NWBUF) CALL GFPART(IPART(1), NAPART, ITRTYP, AMASS, CHARGE, & TLIFE, UBUF, NWBUF) call gfvert(nvert, vert, ntbeam, ngen, tofgn, ubuf, nwbuf) call uhtoc(ibuf(3),4, name1, 4) call uhtoc(ibuf(4),4, name2, 4) c c ubuf is not filled properly for initial vertex.. c if((debug(4).ge.1) .and. (ntbeam .ne.0)) then write(6,*) 'Track = ', indx write(6,*) ibuf(1) write(6,*) ibuf(2) write(6,*) name1 write(6,*) name2 endif if(idebug .and. iswit(4) .eq. ichm) call gdxyz(indx) c Itrack_s(1) = 3 Itrack_s(2) = Tracklist(it)+(abs(nevpar)-1)*10000 Itrack_s(3) = Ipart(1) Itrack_s(4) = charge Rtrack_s(5) = pvert(1) Rtrack_s(6) = pvert(2) Rtrack_s(7) = pvert(3) Rtrack_s(8) = vert(1) Rtrack_s(9) = vert(2) Rtrack_s(10) = vert(3) Itrack_s(11) = ntbeam+(abs(nevpar)-1)*10000 if(ntbeam .gt. 0) then if(nwbuf .gt. 0) then Itrack_s(12) = ibuf(1) else Itrack_s(12) = 0 endif Itrack_s(13) = ibuf(3) Itrack_s(14) = ibuf(4) else Itrack_s(12) = 0 Itrack_s(13) = 0 Itrack_s(14) = 0 endif c c output track if requested c ========================= if(debug(3).eq.1) then write(lun, 1201) Itrack_s(2) 1201 format('** Track No.', i7) write(lun, 1202) Itrack_s(3), Itrack_s(4),Itrack_s(11) 1202 format(' PID :', i8,/, $ ' Charge :', i8,/ $ ' ParentTr :', i8) write(lun, 1203) (Rtrack_s(i), i=5,7) 1203 format(' Momentum :',3f8.2) write(lun, 1204) (Rtrack_s(i), i=8,10) 1204 format(' Vertex :',3f8.2) write(lun, 1205) Itrack_s(12) 1205 format(' ParentVol:', i8) endif c c output this track c if(first_hit) then first_hit = .false. nb = 6*4 nbw = cwrite(fileid, Ievt_s, nb) if(nbw .ne. nb) then write(6,*) 'output error', nb, nbw endif endif nb = 14*4 nbw = cwrite(fileid, Itrack_s, nb) if(nbw .ne. nb) then write(6,*) 'output error', nb, nbw endif end do c c output a -802 if last sub event for this main event... c changed from -1 to -802 to make search for errors easier. c The -1 could be mistaken for a negative chage in the record. c if(nevpar .lt. 0) then iword = -802 nb = 4 nbw = cwrite(fileid, iword, nb) if(nbw .ne. nb) then write(6,*) 'output error', nb, nbw endif endif * 99 RETURN END c c c ========================================= logical Function IsInTrackList(tracknum) c ========================================= Integer tracknum c c see if the 'new' tracknum is already in list c integer num_tracks integer tracklist common /track_list/ num_tracks, tracklist(500) c do i = 1, num_tracks if(tracklist(i) .eq. tracknum) then isintrackList = .true. return endif end do IsInTrackList = .false. return end c c ========================================= subroutine init_detector c ========================================= implicit none c c common block for comunicating c integer max_volume parameter (max_volume=50) character*4 detector_names(max_volume,10) character*4 detector_set(max_volume) integer detector_list(max_volume,20) integer detector_nsub(max_volume) integer number_detectors c common /setup/ $ detector_list, detector_names, $ detector_set, $ detector_nsub, number_detectors c c c local declarations c ================== character*4 sublist(20) character *4 name1, name2 integer ichm, isub, is integer isub_t3, isub_t4,isub_t5 c c Find the level of detector description c integer nlev integer lnam(20), lnum(20), numbv(20) c call saveoption(0) c c Loop over chamber information c ============================= c This is specifically setup for FS. c There need to be flags/ options to include other detectors c e.g. for MRS, Beam-Beam mult c Strategy has to be rethought c ---------------------------- number_detectors = 30 do ichm = 1, number_detectors c isub = 1 IF(ICHM .EQ. 1) THEN NAME1 = 'T1 ' NAME2 = 'T1PR' ELSE IF(ICHM .EQ. 2) THEN NAME1 = 'T2 ' NAME2 = 'T2PR' ELSE IF(ICHM .EQ. 3) THEN NAME1 = 'H1 ' NAME2 = 'TFS1' ELSE IF(ICHM.EQ. 4) THEN isub_t3 = 0 call zfpath('T3 ','T3 ', numbv, nlev, lnam, lnum) if(nlev.eq.0) then call zfpath('T3 ','T3A1', numbv, nlev, lnam, lnum) if(nlev .gt. 0) then isub_t3 = 1 endif endif c c The setup of T3 depends on subvolume setup c The value of isub_t3 has been set in ntinit. c if(isub_t3.eq.0) then NAME1 = 'T3 ' NAME2 = 'T3 ' elseif(isub_t3 .eq. 1) then NAME1 = 'T3 ' NAME2 ='T3XX' isub=3 sublist(1) = 'T3A1' sublist(2) = 'T3A2' sublist(3) = 'T3A3' endif ELSE IF(ICHM .EQ. 5) THEN isub_t4 = 0 call zfpath('T4 ','T4 ', numbv, nlev, lnam, lnum) if(nlev.eq.0) then call zfpath('T4 ','T4A1', numbv, nlev, lnam, lnum) if(nlev .gt. 0) then isub_t4 = 1 endif endif c c The setup of T4 depends on subvolume setup c The value of isub_t4 has been set in ntinit. c if(isub_t4.eq.0) then NAME1 = 'T4 ' NAME2 = 'T4 ' elseif(isub_t4 .eq. 1) then NAME1 = 'T4 ' NAME2 ='T4XX' isub=3 sublist(1) = 'T4A1' sublist(2) = 'T4A2' sublist(3) = 'T4A3' endif ELSE IF(ICHM .EQ. 6) THEN isub_t5 = 0 call zfpath('T5 ','T5 ', numbv, nlev, lnam, lnum) if(nlev.eq.0) then call zfpath('T5 ','T5A1', numbv, nlev, lnam, lnum) if(nlev .gt. 0) then isub_t5 = 1 endif endif c c The setup of T5 depends on subvolume setup c The value of isub_t5 has been set in ntinit. c if(isub_t5.eq.0) then NAME1 = 'T5 ' NAME2 = 'T5 ' elseif(isub_t5 .eq. 1) then NAME1 = 'T5 ' NAME2 ='T5XX' isub=3 sublist(1) = 'T5A1' sublist(2) = 'T5A2' sublist(3) = 'T5A3' endif ELSE IF(ICHM.eq. 7) then NAME1 = 'H2 ' NAME2 = 'TFS2' ELSE IF(ICHM.eq. 8) then NAME1 = 'C1 ' NAME2 = 'C1VA' ELSE IF(ICHM.eq. 9) then NAME1 = 'C3 ' NAME2 = 'RICH' ELSE IF(ICHM .EQ. 10) THEN NAME1 = 'TPM1' NAME2 = 'M1PR' ELSE IF(ICHM .EQ. 11) THEN NAME1 = 'TPM2' NAME2 = 'M2PR' ELSE IF(ICHM .EQ. 12) THEN NAME1 = 'TOFW' NAME2 = 'TFSL' ELSE IF(ICHM .EQ. 13) THEN NAME1 = 'WF1 ' NAME2 = 'WF1 ' ELSE IF(ICHM .EQ. 14) THEN NAME1 = 'WF2 ' NAME2 = 'WF2 ' ELSE IF(ICHM .EQ. 15) THEN NAME1 = 'WM1 ' NAME2 = 'WM1 ' ELSE IF(ICHM .EQ. 16) THEN NAME1 = 'WM2 ' NAME2 = 'WM2 ' c c Beam Beam counters c ================== ELSE IF(ICHM .EQ. 17) THEN NAME1 = 'BEAM' NAME2 = 'PURA' ELSE IF(ICHM .EQ. 18) THEN NAME1 = 'BEAM' NAME2 = 'PURB' ELSEIF(ICHM .EQ. 19) THEN NAME1 = 'BEAM' NAME2 = 'PULA' ELSE IF(ICHM .EQ. 20) THEN NAME1 = 'BEAM' NAME2 = 'PULB' c c New TOFW c ======== ELSE IF(ICHM .EQ. 21) THEN NAME1 = 'TFP1' NAME2 = 'TFT1' else if(ichm .eq. 22) then name1 = 'TILE' name2 = 'TILE' else if(ichm .eq. 23) then name1 = 'STRA' name2 = 'STRA' else if(ichm .eq. 24) then name1 = 'STRB' name2 = 'STRB' else if(ichm .eq. 25) then name1 = 'STRC' name2 = 'STRC' else if(ichm .eq. 26) then name1 = 'STRD' name2 = 'STRD' else if(ichm .eq. 27) then name1 = 'STRE' name2 = 'STRE' else if(ichm .eq. 28) then name1 = 'STRF' name2 = 'STRF' else if(ichm .eq. 29) then name1 = 'STRG' name2 = 'STRG' ELSE IF(ICHM .EQ. 30) THEN NAME1 = 'TFP2' NAME2 = 'TFT2' ENDIF c detector_nsub(ichm) = isub c c Name1 is the set name c Name2 is the physical volume name (active) c if(isub .eq. 1) then sublist(1) = NAME2 endif c c store the subnames c detector_set(ichm) = name1 do is = 1, isub detector_names(ichm,is) = sublist(is) end do end do return end c ========================================= subroutine saveoption(i) c ========================================= integer i logical save_all, save_charged, save_primary common /out_option/ save_all, save_charged, save_primary c c save_all : if(i .eq. 0) then save_all = .true. save_charged = .false. save_primary=.false. elseif(i.eq.1) then save_all = .false. save_charged = .true. save_primary=.false. endif if(i.eq.2) save_primary=.true. return end c ========================================= subroutine gbrinit(fileid) c ========================================= implicit none c integer fileid c integer max_volume parameter (max_volume=50) character*4 detector_names(max_volume,10) character*4 detector_set(max_volume) integer detector_list(max_volume,20) integer detector_nsub(max_volume) integer number_detectors c integer version_id parameter (version_id=1001) c integer isub_t5, isub_t4, isub_t3 common /setup/ $ detector_list, detector_names, $ detector_set, $ detector_nsub, number_detectors c c Find the level of detector description c integer nlev integer lnam(20), lnum(20), numbv(20) c INTEGER NLEVEL,NAMES,NUMBER,LVOLUM,LINDEX,INFROM,NLEVMX, + NLDEV,LINMX REAL GTRAN,GRMAT,GONLY,GLX * COMMON/GCVOLU/NLEVEL,NAMES(15),NUMBER(15), + LVOLUM(15),LINDEX(15),INFROM,NLEVMX,NLDEV(15),LINMX(15), + GTRAN(3,15),GRMAT(10,15),GONLY(15),GLX(3) c character*4 name1, name2 character*4 sublist(20) integer ichm, isub, is, ier, i, j, ivol real param(15) integer npar, ishape, nb, k integer cwrite character *4 magnet_set(5) data magnet_set /'D1 ','D2 ','D3 ','D4 ','D5 '/ integer iname * vector debug(5) * c c This routine will find geometry parameter ans names of active c volumes and will write them to the output file as an Run record c c The run record is initially defined as c c c struct Volume { c int Ishape c int Npar c char name[4]; c float pos[3]; c float RotAng[3][3]; c float Param[12]; c int Active c } c c struct Run { c int NoVolumes; c Volume volumes[NoVolumes] c } c integer volsize parameter (volsize=30) real fVolume(volsize,50) integer IVolume(volsize,50) equivalence (IVolume, fVolume) Integer Nvolumes c integer lun_d c c Some thought has to be put into the dealing of the c beam counter volumes. Since they consist of a many level above c the lower level the GLVOLU etc does not work properly on c the geometry extraction. There is only an active volume c e.g. PURA in PMRA but the hit number comes from PMRA not the lower level c active volume. Is the hit recorded properly ? c c c Loop over chanber information c ============================= NVolumes = 0 do ichm = 1, number_detectors name1 = detector_set(ichm) isub = detector_nsub(ichm) lun_d = 6 do is = 1,isub name2 = detector_names(ichm, is) call zfpath(name1, name2, numbv, nlev, lnam, lnum) if(debug(5).gt.0) then write(lun_d,*) '*** Volume list ***' do k = 1, nlev call uhtoc(lnam(k),4,name1,4) write(lun_d,*) name1, lnum(k) end do endif c in the case of a non existing detector and/or set (nlev <=0) c the remaining part of this loop is skipped c This should though not happen. if(nlev.gt. 0) then if(debug(5).gt.0) then write(6,*) nlev, lnum(nlev) endif if(lnum(nlev).eq.0) nlev=nlev-1 NVolumes = Nvolumes +1 detector_list(ichm,is) = Nvolumes if(name1 .eq. 'TILE') then Ivolume(3,NVolumes) = lnam(nlev+1) else Ivolume(3,NVolumes) = lnam(nlev) endif call glvolu(nlev, Lnam, lnum, ier) call uhtoc(lnam(nlev),4,name2,4) call GFVOLU(name2, ivol, ishape, npar, param) Ivolume(1,NVolumes) = ishape Ivolume(2,NVolumes) = npar do k=1,npar FVolume(15+k,NVolumes) = param(k) end do do k= npar+1, 12 FVolume(15+k, NVolumes) = 0.0 end do IVolume(28, NVolumes) = 1 do i=0,2 do k=0,2 FVolume(7+k+i*3, NVolumes) = grmat(k*3+i+1, nlev) end do end do do i=1,3 FVolume(3+i,NVolumes) = GTRAN(i,nlev) end do if(debug(5).gt.0) then write(6,*) 'Volume number : ', ivol write(6,*) 'Shape :',ishape, npar write(6, 99) (param(i),i=1,npar) 99 format(1x,3(1x,f7.3)) write(6,*) 'Set : ', name1,' Volume : ', name2 write(6,*) 'Position' write(6,100) (GTRAN(i,nlev),i=1,3) 100 format(3(1x,f8.3)) write(6,*) 'Rotation matrix' do j=1,3 write(6,101) (grmat((j-1)*3+i,nlev),i=1,3) 101 format(3(2x,f8.5)) end do endif endif end do end do do i=1,6 lnum(i) = 1 end do c c Set information for magnetic volumes c do is = 1, 5 name1 = magnet_set(is) if(name1 .eq. 'D1 ') then call uctoh('CAVE',lnam(1), 4, 4) call uctoh('FMS1',lnam(2), 4, 4) call uctoh('D1 ',lnam(3), 4, 4) call uctoh('MAG1',lnam(4), 4, 4) call uctoh('GAP1',lnam(5), 4, 4) nlev = 5 elseif(name1 .eq. 'D2 ') then call uctoh('CAVE',lnam(1), 4, 4) call uctoh('FMS1',lnam(2), 4, 4) call uctoh('FD2 ',lnam(3), 4, 4) call uctoh('D2 ',lnam(4), 4, 4) call uctoh('MAG2',lnam(5), 4, 4) call uctoh('GAP2',lnam(6), 4, 4) nlev = 6 elseif(name1 .eq. 'D3 ') then call uctoh('CAVE',lnam(1), 4, 4) call uctoh('FMS2',lnam(2), 4, 4) call uctoh('D3 ',lnam(3), 4, 4) call uctoh('MAG3',lnam(4), 4, 4) call uctoh('GAP3',lnam(5), 4, 4) nlev = 5 elseif(name1 .eq. 'D4 ') then call uctoh('CAVE',lnam(1), 4, 4) call uctoh('FMS2',lnam(2), 4, 4) call uctoh('D4 ',lnam(3), 4, 4) call uctoh('MAG4',lnam(4), 4, 4) call uctoh('GAP4',lnam(5), 4, 4) nlev = 5 elseif(name1 .eq. 'D5 ') then call uctoh('CAVE',lnam(1), 4, 4) call uctoh('MIDS',lnam(2), 4, 4) call uctoh('M0 ',lnam(3), 4, 4) call uctoh('MAGN',lnam(4), 4, 4) call uctoh('AGAP',lnam(5), 4, 4) nlev = 5 endif if(debug(5).gt.1)then write(lun_d,*) '*** Volume list ***' do k = 1, nlev call uhtoc(lnam(k),4,name1,4) write(lun_d,*) name1, lnum(k) end do endif call glvolu(nlev, lnam, lnum,ier) if(ier .eq. 0) then call uctoh(name1,iname,4,4) lun_d = 6 NVolumes = Nvolumes +1 Ivolume(3,NVolumes) = iname call GFVOLU(name1, ivol, ishape, npar, param) Ivolume(1,NVolumes) = ishape Ivolume(2,NVolumes) = npar do k=1,npar FVolume(15+k,NVolumes) = param(k) end do do k= npar+1, 12 FVolume(15+k,NVolumes) = 0.0 end do IVolume(28,NVolumes) = 0 do i=0,2 do k=0,2 FVolume(7+k+i*3,NVolumes) = grmat(k*3+i, nlev) end do end do do i=1,3 FVolume(3+i,NVolumes) = GTRAN(i,nlev) end do if(debug(5).gt.1) then write(6,*) 'Volume number : ', ivol write(6,*) 'Shape :',ishape, npar write(6, 99) (param(i),i=1,npar) write(6,*) 'Magnet : ', name1 write(6,*) 'Position' write(6,100) (GTRAN(i,nlev),i=1,3) write(6,*) 'Rotation matrix' do j=1,3 write(6,101) (grmat((j-1)*3+i,nlev),i=1,3) end do endif else write(6,*) 'Magnet Volume not found' do k = 1, nlev call uhtoc(lnam(k),4,name1,4) write(lun_d,*) name1, lnum(k) end do endif enddo c c Write the data out c if(debug(5).gt.0) then write(6,*) 'N volumes', NVolumes endif * nb = cwrite(fileid, version_id,4) nb = cwrite(fileid, NVolumes,4) do i=1, NVolumes nb = cwrite(fileid, IVolume(1,i), 112) end do return end SUBROUTINE gbrfile(file) integer o_flag, o_mode character *(*) file * integer copen, cmodeopen * integer fileid common /c_out/ fileid * * open for write+create * tested and working from geant comis call on 12/31/96 * o_flag = cmodeopen() o_mode = 436 l = lenocc(file) fileid = copen(file(1:l), o_flag, o_mode) c ************************************************************************ write(6,*) 'Opening file ',file(1:l), fileid call gbrinit(fileid) 99 RETURN END * SUBROUTINE gbrend * ************************************************************************ integer fileid common /c_out/ fileid if(fileid .ne.0) then write(6,*) 'closing file', fileid call cclose(fileid) endif 99 RETURN END * * subroutine reset_nlevel INTEGER NLEVEL,NAMES,NUMBER,LVOLUM,LINDEX,INFROM,NLEVMX, + NLDEV,LINMX REAL GTRAN,GRMAT,GONLY,GLX * COMMON/GCVOLU/NLEVEL,NAMES(15),NUMBER(15), + LVOLUM(15),LINDEX(15),INFROM,NLEVMX,NLDEV(15),LINMX(15), + GTRAN(3,15),GRMAT(10,15),GONLY(15),GLX(3) nlevel = 0 return end * * Modified versino of GFVOLU (meant to get hold of the * volume parameters. * Subroutine GFVOLU(IUDET,IVOLU, ISHAPE, NPAR, PARAM) Implicit None c c Global Specifications:- c ======================= INTEGER NMATE ,NVOLUM,NROTM,NTMED,NTMULT,NTRACK,NPART + ,NSTMAX,NVERTX,NHEAD,NBIT ,NALIVE,NTMSTO COMMON/GCNUM/NMATE ,NVOLUM,NROTM,NTMED,NTMULT,NTRACK,NPART + ,NSTMAX,NVERTX,NHEAD,NBIT COMMON /GCNUMX/ NALIVE,NTMSTO INTEGER IQ,LQ,NZEBRA,IXSTOR,IXDIV,IXCONS,LMAIN,LR1 INTEGER KWBANK,KWWORK,IWS REAL GVERSN,ZVERSN,FENDQ,WS,Q C PARAMETER (KWBANK=69000,KWWORK=5200) COMMON/GCBANK/NZEBRA,GVERSN,ZVERSN,IXSTOR,IXDIV,IXCONS,FENDQ(16) + ,LMAIN,LR1,WS(KWBANK) integer big parameter (big=16000000) DIMENSION IQ(big),Q(big),LQ(big),IWS(2) EQUIVALENCE (Q(1),IQ(1),LQ(9)),(LQ(1),LMAIN),(IWS(1),WS(1)) INTEGER JDIGI ,JDRAW ,JHEAD ,JHITS ,JKINE ,JMATE ,JPART + ,JROTM ,JRUNG ,JSET ,JSTAK ,JGSTAT,JTMED ,JTRACK,JVERTX + ,JVOLUM,JXYZ ,JGPAR ,JGPAR2,JSKLT C COMMON/GCLINK/JDIGI ,JDRAW ,JHEAD ,JHITS ,JKINE ,JMATE ,JPART + ,JROTM ,JRUNG ,JSET ,JSTAK ,JGSTAT,JTMED ,JTRACK,JVERTX + ,JVOLUM,JXYZ ,JGPAR ,JGPAR2,JSKLT C c c Sub-program argument specifications:- c ===================================== Character*4 IUDET ! Detector name Integer IVOLU Integer ISHAPE, NPAR real PARAM(*) c c local declarations c ================== Integer JVO, NIN, JVOMOT,MNIN,NMBR,I,KONLY,MOTHER,IVOMGL,JIN,NATT $ ,ATT,INGLOB,jpar character *4 name c c Executable code c c =============== CALL GLOOK(IUDET,IQ(JVOLUM+1),NVOLUM,IVOLU) JVO = LQ(JVOLUM-IVOLU) ISHAPE = Q(JVO+2) NIN = Q(JVO+3) NPAR = Q(JVO+5) NATT = Q(JVo+6) JPAR = JVO + 6 CALL UCOPY (Q(JPAR+1), PARAM, NPAR) return end c =========================================================== Subroutine ZFPATH (IUSET, IUDET, NUMBV, NLEV, LNAME, LNUMB) c =========================================================== c c Description:- c ============= c c Routine to augment GFPATH. Called with detector and set c names not indices. GFPATH returns the full path through c to volume tree for sensitive detectors defined using c GSDETV. NUMBV should contain only indeterminate c volume numbers (i.e. for levels at which there are c multiple copies of a volume). LNAME and LNUMB contain c the full path through the volume tree on return, and c NLEV contains the number of levels in the full path. c c Note: If ZFPATH fails to find the SET bank, a value of c ----- of NLEV=-1 is returned. If ZFPATH fails to find c the the DET bank, NLEV=0 is returned. c c Author: BAC c =========== c c Creation Date: 22-SEP-1992 c ========================== c c Modification history c ==================== c c Implicit inputs, outputs, side effects:- c ======================================== Implicit None c c Global Specifications:- c ======================= INTEGER IQ,LQ,NZEBRA,IXSTOR,IXDIV,IXCONS,LMAIN,LR1 INTEGER KWBANK,KWWORK,IWS REAL GVERSN,ZVERSN,FENDQ,WS,Q C PARAMETER (KWBANK=69000,KWWORK=5200) COMMON/GCBANK/NZEBRA,GVERSN,ZVERSN,IXSTOR,IXDIV,IXCONS,FENDQ(16) + ,LMAIN,LR1,WS(KWBANK) integer big parameter (big=16000000) DIMENSION IQ(big),Q(big),LQ(big),IWS(2) EQUIVALENCE (Q(1),IQ(1),LQ(9)),(LQ(1),LMAIN),(IWS(1),WS(1)) INTEGER JDIGI ,JDRAW ,JHEAD ,JHITS ,JKINE ,JMATE ,JPART + ,JROTM ,JRUNG ,JSET ,JSTAK ,JGSTAT,JTMED ,JTRACK,JVERTX + ,JVOLUM,JXYZ ,JGPAR ,JGPAR2,JSKLT C COMMON/GCLINK/JDIGI ,JDRAW ,JHEAD ,JHITS ,JKINE ,JMATE ,JPART + ,JROTM ,JRUNG ,JSET ,JSTAK ,JGSTAT,JTMED ,JTRACK,JVERTX + ,JVOLUM,JXYZ ,JGPAR ,JGPAR2,JSKLT C c c Sub-program argument specifications:- c ===================================== Character*4 IUSET ! Set name Character*4 IUDET ! Detector name Integer NUMBV (*) ! Detector volume numbers Integer NLEV ! Number of volume tree levels for det. Integer LNAME (*) ! Volume tree names for det. Integer LNUMB (*) ! Volume tree numbers for det. c c External Specifications:- c ========================= c c Local Specifications:- c ====================== Integer JS , ! Link into detector banks for set & NDET , ! Number of defined hit sets & NSET , ! Number of detectors in set & IDET , ! Detector index & ISET ! Set index c c Executable Statements:- c ======================= c c Find Set index from list of GEANT set names c =========================================== NSET = IQ (JSET - 1) CALL GLOOK (IUSET, IQ (JSET + 1), NSET, ISET) IF (ISET .LE. 0) Then NLEV = -1 Go to 999 End if c c Find link to set bank c ===================== JS = LQ (JSET - ISET) IF (JS .LE. 0) Then NLEV = 0 Go to 999 End if c c Find detector index from list of GEANT det. names c ================================================= NDET = IQ (JS - 1) CALL GLOOK (IUDET, IQ (JS + 1), NDET, IDET) IF (IDET .LE. 0) Then NLEV = 0 Go to 999 End if c c Now call GFPATH using found set and det. indices c ================================================ Call GFPATH (ISET, IDET, NUMBV, NLEV, LNAME, LNUMB) 999 Continue Return End c ============================================ subroutine print_hits(lun, setname, detname, $ ihits_s, rhits_s) c ============================================ c integer lun integer ihits_s(21) character*4 setname character*4 detname c real rhits_s(21) c output hit if requested c ========================= write(lun, 1101) setname, detname 1101 format('** HIT from set ',a4,' detector ' , a4) write(lun,1102) ihits_s(2), ihits_s(3), ihits_s(20) 1102 format(' Track no: ', I8,/, $ ' Detector: ', i8,/, $ ' Pid: ',i8) write(lun,1103) (rhits_s(i),i=4,6) 1103 format(1x,'LocalIn :',3f8.2) write(lun,1104) (rhits_s(i),i=7,9) 1104 format(1x,'LocalOut :',3f8.2) write(lun,1105) (rhits_s(i),i=10,12) 1105 format(1x,'GlobalIn :',3f8.2) write(lun,1106) (rhits_s(i),i=13,15) 1106 format(1x,'GlobalOut:',3f8.2) write(lun,1107) rhits_s(16)*1.0E9, $ rhits_s(17)*1.0E3, $ rhits_s(18) 1107 format(1x,'Tof :',f8.2,/ $ ' de/dx :',f8.2,/ $ ' pDet :',f8.2) write(lun,1108) ihits_s(19) 1108 format(1x,'SubVolNo :',i8) return end