* * $Id: tpc_hitlist.F,v 1.6 2004/05/28 16:06:39 dpp Exp $ * * $Log: tpc_hitlist.F,v $ * Revision 1.6 2004/05/28 16:06:39 dpp * -> use input codes * -> fix error is recognizing valid hits at bottom of curler * -> fix error in assigning layer namber at bottom and top of curler * -> allow fragmented hit list * -> skip processing hits for L_EVENT_MATCH=F * -> skip printing hits for KILL=1 * -> dialog to input overide values for windo limits * -> add variable for cdscrtcd hit usage * * Revision 1.5 2004/01/23 21:11:29 dpp * -> use LCD crossings as borders of cells * -> create multiple ionization centers for input into * tpc_pad_response * -> add MC tag to argument of tpc_ionization centers * -> initialize the tag match value * -> initialize event match * -> add logical argument to indicate event match * * Revision 1.4 2003/10/28 21:26:07 dpp * -> diagnostics * -> rename doit variable "hittrk" to HIT2TRK * -> rename doit variable "trkhit" to TRK2HIT * -> identify list of plausible mc tracks * -> use tpc_hitlist_xcheck to find number of isolated hits on a * plausible mc track * * Revision 1.3 2003/08/28 13:07:09 dpp * -> pointer to first instance of a device in the event hit list * * Revision 1.2 2003/05/12 16:08:27 dpp * -> use CLEO_DOITTPC to dummy this routine in CLEO production * * Revision 1.1 2003/03/05 17:50:35 dpp * -> NEW ROUTINE * -> main routine to convert LCD hit list to cdgeom * -> changed all variable names in cdscrtcd to have common root * * #include "sys/CLEO_machine.h" #include "pilot.h" SUBROUTINE TPC_HITLIST(L_EVENT_MATCH) C....................................................................... C. C. TPC_HITLIST - transfer hit list from LCD into cdgeom C. C. COMMON : C. CALLS : C. CALLED : C. AUTHOR : D. Peterson C. C. VERSION : 1.00 C. CREATED : 22-Oct-2002 C. C....................................................................... #if defined(CLEO_TYPCHK) IMPLICIT NONE #endif SAVE LOGICAL L_EVENT_MATCH #if defined(CLEO_DOITTPC) #include "doit/duseq/controlpar.inc" #include "doit/duseq/tfconspa.inc" #include "doit/duseq/tfindpar.inc" C which includes C #include "cl3seq/cdgm3/cdgeompa.inc" C #include "doit/duseq/tfgeompa.inc" #include "doit/duseq/tfctlcde.inc" #include "doit/duseq/tfgeomcd.inc" #include "cl3seq/cdgm3/cdgeomcd.inc" #include "doit/duseq/cdscrtcd.inc" #include "doit/duseq/runev.inc" #include "doit/sfseq/sfpar.inc" #include "doit/sfseq/sfcom.inc" #include "doit/sfseq/sfxtsp.inc" #include "doit/sfseq/sfxtsc.inc" #include "doit/duseq/tftrakcd.inc" #include "doit/duseq/monte_carlo_tracks.inc" #include "doit/duseq/tpccom.inc" #include "/home/dpp/lcd_simulation/cornell/hep/lcd/io/fortran/lcdevt.inc" INTEGER print_EvntHitMax,print_limitofcells,print_limitoflayers INTEGER I,J INTEGER IGO,KILL INTEGER TRK_NOW,TRK_OLD LOGICAL VALID_TPC_CROSSING INTEGER ICROSS REAL XCROSS,YCROSS INTEGER ICROSS_ENTR,ICROSS_2PREV,ICROSS_ENTR_H,ICROSS_2PREV_H REAL RAD_EXIT,PHI_EXIT,Z_EXIT REAL RAD_ENTR,PHI_ENTR,Z_ENTR LOGICAL CREATE_HITS_WITH_CROSSING REAL ATN2PI REAL Z_NOW REAL Z_STRT,Z_STEP INTEGER ILCD INTEGER IFIT REAL RADSQ REAL XTMP,YTMP,FTMP INTEGER IERROR REAL ARCFUN REAL CALCPPERP,CALCMOM(3) INTEGER FIRST_LYR_SAVED INTEGER LAST_LYR_SAVED INTEGER FIRST_HIT_SAVED INTEGER LAST_HIT_SAVED LOGICAL ACCEPT_MORE_HITS INTEGER NCLEAR INTEGER MC_temp,NH_temp,NC_temp,it,jt,imax,jmin REAL CU_temp,FI_temp,D0_temp,DZ_temp,Z0_temp LOGICAL FIRST_ENTRY/.TRUE./ c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 * ----------Executable code starts here--------------------------------- C----------------------------------------------------------------------- C event information C----------------------------------------------------------------------- current_event=current_event+1 L_EVENT_MATCH=.TRUE. print 1000,current_event,NTRKHITS,NMCPART,NEVLCD 1000 FORMAT( 1 ' TPC_HITLIST: seq event:',I8,' found',I9,' crossings'/ 2 32X,' found',I9,' MC particles'/ 3 19X,' found event number',I9) C----------------------------------------------------------------------- #if defined(CLEO_SFDIAG) IF(FIRST_ENTRY)THEN KILL=0 TAG_MATCH=0 EVENT_MATCH=0 FIRST_ENTRY=.FALSE. print 1002 1002 FORMAT(' TPC_HITLIST: enter the EVENT NUMBER to display') read *,EVENT_MATCH print 1003,EVENT_MATCH 1003 format(' will stop at EVENT',I8) print 1004 1004 FORMAT(' TPC_HITLIST: enter the track ID to tag in display') read *,TAG_MATCH print 1005,TAG_MATCH 1005 format(' will show track ID',I8,' in different color') print 1016 1016 format(' enter lower window limit; -999 will inhibit overide') read *,SFZSLWLOVRD if(SFZSLWLOVRD.le.-999)then SFZSLWLOVRD=-999 SFZSLWMOVRD=999 else if(SFZSLWLOVRD.lt.-80)SFZSLWLOVRD=-80 if(SFZSLWLOVRD.gt.80)SFZSLWLOVRD=80 print 1017 1017 format(' enter UPPER window limit; will force >/= lower') read *,SFZSLWMOVRD if(SFZSLWMOVRD.gt.80)SFZSLWMOVRD=80 if(SFZSLWMOVRD.lt.SFZSLWLOVRD)SFZSLWMOVRD=SFZSLWLOVRD endif ENDIF L_EVENT_MATCH=(current_event.EQ.EVENT_MATCH) IF(L_EVENT_MATCH)THEN KILL=0 ENDIF #else KILL=0 TAG_MATCH=1 EVENT_MATCH=1 SFZSLWLOVRD=-999 SFZSLWMOVRD=999 #endif print_EvntHitMax=EvntHitMax EvntHit_Num=0 print 1011,print_EvntHitMax 1011 FORMAT(' TPC_HITLIST: limit of number of hits=',I9, 2 ' zero number of hits: EvntHit_Num') C initilize the hit storage location link list EvntHit_OpnLoc1=1 EvntHit_OpnLocS=1 EvntHit_OpnLocL=EvntHitMax DO 53 I=EvntHit_OpnLoc1,EvntHit_OpnLocL EvntHit_OpnLocN(I)=I+1 EvntHit_OpnLocP(I)=I-1 53 CONTINUE EvntHit_OpnLocN(EvntHit_OpnLocL)=0 EvntHit_OpnLocP(EvntHit_OpnLoc1)=0 print 1012 1012 FORMAT(' TPC_HITLIST: initialize the 1st open location;', 2 ' initialize the next open locations') print_limitofcells=KWS1CD CALL VZERO(EvntHit_MapDet,KWS1CD) print 1013,print_limitofcells 1013 FORMAT(' TPC_HITLIST: limit of number of cells=',I9, 2 ' zero the hit map: EvntHit_MapDet') print_limitoflayers=KLYRCD CALL VZERO(EvntHit_1stLyr,KLYRCD) print 1014,print_limitoflayers 1014 FORMAT(' TPC_HITLIST: limit of number of layers=',I9, 2 ' zero the 1st hit pointer: EvntHit_1stLyr') CALL VZERO(EvntHit_1stDEV,MXDVCD) n_mctrak_plausbl=0 C----------------------------------------------------------------------- C process the input hit list C loop over hits C----------------------------------------------------------------------- trk_old=-2 IF( 1 (L_EVENT_MATCH).AND. 2 (NTRKHITS.GT.1))THEN EvntHit_1stDEV(ITPC1)=1 do 219 i=1,NTRKHITS C----------------------------------------------------------------------- C recognize new track C----------------------------------------------------------------------- C There is a miss-match beteen the particle momentum and the hits C that I do not understand but has something to do with Java counting C from 0 instead of 1. So, for now, put in a cludge to the C MC particle number. TRK_NOW=trkhitmcpart(i) +1 IF(TRK_NOW.ne.trk_old)THEN C----------------------------------------------------------------------- C finish up previous track C possibly store in list of plausible MC tracks, those that could be found C----------------------------------------------------------------------- IF( 1 (trk_old.GT.0).and. 2 (NFIT.GT.6).and. 3 (FIRST_LYR_SAVED.EQ.1).and. 4 (LAST_LYR_SAVED.GE.30).and. 5 (n_mctrak_plausbl.LT.M_mctrak_plausbl) 6 )THEN print 2000,trk_old,NFIT,first_lyr_saved,last_lyr_saved 2000 format(' tpc_hitlist:will fit input track:', 1 ' number',I5,' NFIT=',I3, 2 ' first layer saved=',I4, 2 ' last layer saved=',I4) CALL CFCFIT PHI6CF=0. RKN6CF=999. print 2011,trk_old,NFIT,KAPCF,ALPCF,BETCF,GAMCF,XICF 2011 format(' tpc_hitlist: fitted input track:', 1 ' number ',I5,' NFIT=',I3, 2 ' K=',F10.5,' A=',F10.5,' B=',F10.5, 3 ' G=',F10.5,' XI=',F10.5) CALL CDCFID( 1 CURCF,PHI0CF,D0CF, 2 KAPCF,ALPCF,BETCF,GAMCF,XICF) print 2012,CURCF,PHI0CF,D0CF 2012 format(' tpc_hitlist: fitted input track:',22X, 2 ' C=',F10.5,' F=',F10.5,' D=',F10.5) c must get the arc length DO 119 IFIT=1,NFIT RADSQ=XFIT(IFIT)**2+YFIT(IFIT)**2 CALL TFPHTR(ALPCF,BETCF,KAPCF,GAMCF,PHI6CF,RKN6CF, 2 RADSQ,+1,XTMP,YTMP,FTMP,IERROR) C NEED TO CHECK IERROR SFIT(IFIT)=ARCFUN(KAPCF,ALPCF,BETCF,XTMP,YTMP,XICF) 119 CONTINUE Z0BIAS=9.9 CALL SZ_FITTER print 2013,trk_old,NFIT,Y0LF,TANDLF 2013 format(' tpc_hitlist: fitted input track:', 1 ' number ',I5,' NFIT=',I3, 2 ' Z0=',F10.5,' DZ=',F10.5) print 2014,trk_old, 1 MCPID(trk_old), 2 MCPMOM(1,trk_old), 3 MCPMOM(2,trk_old), 4 MCPMOM(3,trk_old) 2014 format( 1 ' recall, this was MC particle numb ',I7, 2 ' ID=',I12, 3 ' 3mom=',3F11.3) CALCPPERP=3.*0.15/CURCF CALCMOM(1)=CALCPPERP*COS(PHI0CF) CALCMOM(2)=CALCPPERP*SIN(PHI0CF) CALCMOM(3)=CALCPPERP*TANDLF print 2015, 2 CALCMOM(1), 3 CALCMOM(2), 4 CALCMOM(3) 2015 format( 1 ' calculated momentum from fitted crossings', 2 ' ',15X 3 ' 3mom=',3F11.3) n_mctrak_plausbl=n_mctrak_plausbl+1 MC_mctrak_plausbl(n_mctrak_plausbl)=trk_old CU_mctrak_plausbl(n_mctrak_plausbl)=CURCF FI_mctrak_plausbl(n_mctrak_plausbl)=PHI0CF D0_mctrak_plausbl(n_mctrak_plausbl)=D0CF DZ_mctrak_plausbl(n_mctrak_plausbl)=TANDLF Z0_mctrak_plausbl(n_mctrak_plausbl)=Y0LF NH_mctrak_plausbl(n_mctrak_plausbl)= 2 LAST_HIT_SAVED-FIRST_HIT_SAVED+1 CALL TPC_HITLIST_XCHECK(trk_old, 2 FIRST_HIT_SAVED,LAST_HIT_SAVED, 3 NCLEAR) NC_mctrak_plausbl(n_mctrak_plausbl)=NCLEAR ENDIF C----------------------------------------------------------------------- C ^- end of processing finish of previous track C----------------------------------------------------------------------- C----------------------------------------------------------------------- C decide to READ IN the new track C----------------------------------------------------------------------- print 3013,TRK_NOW, 1 MCPID(TRK_NOW), 2 MCPMOM(1,TRK_NOW), 3 MCPMOM(2,TRK_NOW), 4 MCPMOM(3,TRK_NOW) 3013 format(' ~~~~~~~~~~~'/ 1 ' recognized new MC particle number ',I7, 2 ' ID=',I12, 3 ' 3mom=',3F11.3) #if defined(CLEO_SFDIAG) IF(kill.eq.0)THEN print 3014,TRK_NOW 3014 format( 1 ' for hits on track',I7, 2 ' "N" to jump,', 3 ' or "K" to stop asking') call dsf_rfv_input(igo) if(igo.eq.DSFRFVINPUT_No )go to 220 if(igo.eq.DSFRFVINPUT_Kill)kill=1 ENDIF #endif C----------------------------------------------------------------------- C initialization for new track C----------------------------------------------------------------------- C----------------------------------------------------------------------- C initialization for hits_on_pads generation C----------------------------------------------------------------------- ICROSS_ENTR=-1 ICROSS_2PREV=-1 C----------------------------------------------------------------------- C initialization for input track fit and acceptance for expected find C----------------------------------------------------------------------- NFIT=0 FIRST_LYR_SAVED=0 LAST_LYR_SAVED=0 ACCEPT_MORE_HITS=.TRUE. C----------------------------------------------------------------------- C end of condition: recognize new track number C----------------------------------------------------------------------- trk_old=TRK_NOW ENDIF C----------------------------------------------------------------------- C process the input crossing C----------------------------------------------------------------------- ICROSS=trkhitlayer(i) XCROSS=trkhit(1,i)/100. YCROSS=trkhit(2,i)/100. RAD_EXIT=SQRT(XCROSS**2+YCROSS**2) PHI_EXIT=ATN2PI(YCROSS,XCROSS) Z_EXIT=trkhit(3,i)/100. C----------------------------------------------------------------------- C hold these for the print, below C----------------------------------------------------------------------- ICROSS_ENTR_H=ICROSS_ENTR ICROSS_2PREV_H=ICROSS_2PREV C----------------------------------------------------------------------- C require match of hit radius to layer radius C this eliminates the intermediate detector hits C C When using the LCD crossings as TPC layers boundries, C trkhitlayer is numbered 0:NLYRCD C Thus, trkhitlayer(i) equals ILCD of the TPC layer below it. C When trkhitlayer(i) equals 0, this is the crossing inside the first layer. C----------------------------------------------------------------------- VALID_TPC_CROSSING=.FALSE. CREATE_HITS_WITH_CROSSING=.FALSE. IF(TRK_NOW.GT.0)THEN IF(ICROSS.GT.0)THEN IF(abs(RAD_EXIT-(RCD(ICROSS)+CELRCD(1,ICROSS)/2.)) 1 .lt.0.005)VALID_TPC_CROSSING=.TRUE. ELSEIF(ICROSS.EQ.0)THEN IF(abs(RAD_EXIT-(RCD(1)-CELRCD(1,1)/2.)) 1 .lt.0.005)VALID_TPC_CROSSING=.TRUE. ENDIF ENDIF C----------------------------------------------------------------------- C processing of valid crossings C----------------------------------------------------------------------- IF(VALID_TPC_CROSSING)THEN C can't do anything without a valid previous crossing IF(ICROSS_ENTR.GE.0)THEN C outgoing IF(ICROSS.eq.(ICROSS_ENTR+1))THEN CREATE_HITS_WITH_CROSSING=.TRUE. ILCD=ICROSS C incoming ELSEIF(ICROSS.eq.(ICROSS_ENTR-1))THEN CREATE_HITS_WITH_CROSSING=.TRUE. ILCD=ICROSS+1 C crossing within a layer C sorry, can't process an initial double crossing in below layer 1 ELSEIF( 1 (ICROSS.eq.ICROSS_ENTR).and. 2 (ICROSS_2PREV.GE.0) 3 )THEN C top of curler, ignore escape from chamber IF(ICROSS.EQ.(ICROSS_2PREV+1))THEN IF(ICROSS.LT.NLYRCD)THEN CREATE_HITS_WITH_CROSSING=.TRUE. ILCD=ICROSS+1 ENDIF C bottom of curler, ignore escape from chamber ELSEIF(ICROSS.EQ.(ICROSS_2PREV-1))THEN IF(ICROSS.GT.0)THEN CREATE_HITS_WITH_CROSSING=.TRUE. ILCD=ICROSS ENDIF ENDIF ENDIF ENDIF C----------------------------------------------------------------------- C create hits when there are valid entry and exit crossings C----------------------------------------------------------------------- IF(CREATE_HITS_WITH_CROSSING)THEN IF(ILCD.LE.ARTIFICIAL_LAYER_STOP)THEN c print 9901,ilcd, c 2 PHI_EXIT,PHI_ENTR,Z_EXIT,Z_ENTR,TRK_NOW c 9901 format(' tpc_hitlist: diag:',' ilcd=',I3, c 1 ' phi_exit=',F7.3,' phi_enrt=',F7.3, c 2 ' Z_exit=',F7.3,' Z_enrt=',F7.3, c 3 ' mcpart=',I6) CALL TPC_IONIZATION_CENTERS(ILCD, 2 PHI_EXIT,PHI_ENTR,Z_EXIT,Z_ENTR,TRK_NOW) ENDIF C----------------------------------------------------------------------- C fill arrays to do a fit of the hits, obtain the generated track parameters C----------------------------------------------------------------------- IF(ACCEPT_MORE_HITS)THEN IF(FIRST_LYR_SAVED.EQ.0)THEN FIRST_LYR_SAVED=ILCD FIRST_HIT_SAVED=I ELSE IF(ILCD.LE.LAST_LYR_SAVED)ACCEPT_MORE_HITS=.FALSE. ENDIF ENDIF IF(ACCEPT_MORE_HITS)THEN LAST_LYR_SAVED=ILCD LAST_HIT_SAVED=I C there is room for KLR2TF hits in the track arrays C IADFIT must be faked C CFCFIT does not use IADFIT C CFMTRX requires non-zero IADFIT C SZ_FITTER requires non-zero IADFIT C HIT2TRK and TRK2HIT not used in the fitters NFIT=NFIT+1 IADFIT(NFIT)= i IPLFIT(NFIT)= ILCD XFIT (NFIT)= trkhit(1,i)/100. YFIT (NFIT)= trkhit(2,i)/100. SFIT (NFIT)= RAD_EXIT ZFIT (NFIT)= trkhit(3,i)/100. DFIT (NFIT)= 0. MESFIT(NFIT)= 0. STRFFT(NFIT)= 2 STSZFT(NFIT)= 2 WGTFIT(NFIT)= 1./(.0001)**2 ENDIF C----------------------------------------------------------------------- C end of condition CREATE_HITS_WITH_CROSSING C----------------------------------------------------------------------- ENDIF C----------------------------------------------------------------------- C end of processing input crossing C save "exit" as "entering" for use with the next input crossing C "2prev" is used only to resolve the ambiguous direction of a track C that starts at the top/bottom of a loop C----------------------------------------------------------------------- ICROSS_2PREV=ICROSS_ENTR ICROSS_ENTR=ICROSS Z_ENTR=Z_EXIT RAD_ENTR=RAD_EXIT PHI_ENTR=PHI_EXIT C----------------------------------------------------------------------- C end of condition VALID_TPC_CROSSING C----------------------------------------------------------------------- ELSE ICROSS_2PREV=-1 ICROSS_ENTR=-1 ENDIF C----------------------------------------------------------------------- C print status, including VALID_TPC_CROSSING, CREATE_HITS_WITH_CROSSING C----------------------------------------------------------------------- IF(KILL.NE.1)THEN c print 3015,i,EvntHit_Num, c 1 ICROSS,ICROSS_ENTR_H,ICROSS_2PREV_H, c 2 VALID_TPC_CROSSING,CREATE_HITS_WITH_CROSSING, c 3 TRK_NOW,(trkhit(j,i),j=1,3), c 4 RAD_EXIT,PHI_EXIT,Z_EXIT 3015 format(' TPC_HITLIST:',I7,I8, 1 ' lyr=',I3,' prev=',I3,' 2prev=',I3, 2 ' v',L1,' c',L1, 3 ' mcpart=',I6,' xyz=',3F7.2, 4 ' RfZ=',F7.4,F6.3,F9.5) ENDIF C----------------------------------------------------------------------- C end of loop over hits C----------------------------------------------------------------------- 219 continue C----------------------------------------------------------------------- C jump point to stop processing of the tracks C----------------------------------------------------------------------- 220 CONTINUE TPC_NhitsGen_Cross=EvntHit_Num CALL TPC_RANDOM_NOISE TPC_NhitsGen_wNoise=EvntHit_Num CALL TPC_FADC_RESPONSE(0,0) TPC_NhitsGen_aFADC=EvntHit_Num C----------------------------------------------------------------------- C reorder the found tracks C----------------------------------------------------------------------- if(n_mctrak_plausbl.gt.1)then imax=n_mctrak_plausbl-1 do 316 it=1,imax jmin=it+1 do 315 jt=jmin,n_mctrak_plausbl if(DZ_mctrak_plausbl(jt).lt.DZ_mctrak_plausbl(it))then MC_temp=MC_mctrak_plausbl(it) NH_temp=NH_mctrak_plausbl(it) NC_temp=NC_mctrak_plausbl(it) CU_temp=CU_mctrak_plausbl(it) FI_temp=FI_mctrak_plausbl(it) D0_temp=D0_mctrak_plausbl(it) DZ_temp=DZ_mctrak_plausbl(it) Z0_temp=Z0_mctrak_plausbl(it) MC_mctrak_plausbl(it)=MC_mctrak_plausbl(jt) NH_mctrak_plausbl(it)=NH_mctrak_plausbl(jt) NC_mctrak_plausbl(it)=NC_mctrak_plausbl(jt) CU_mctrak_plausbl(it)=CU_mctrak_plausbl(jt) FI_mctrak_plausbl(it)=FI_mctrak_plausbl(jt) D0_mctrak_plausbl(it)=D0_mctrak_plausbl(jt) DZ_mctrak_plausbl(it)=DZ_mctrak_plausbl(jt) Z0_mctrak_plausbl(it)=Z0_mctrak_plausbl(jt) MC_mctrak_plausbl(jt)=MC_temp NH_mctrak_plausbl(jt)=NH_temp NC_mctrak_plausbl(jt)=NC_temp CU_mctrak_plausbl(jt)=CU_temp FI_mctrak_plausbl(jt)=FI_temp D0_mctrak_plausbl(jt)=D0_temp DZ_mctrak_plausbl(jt)=DZ_temp Z0_mctrak_plausbl(jt)=Z0_temp endif 315 continue 316 continue endif C----------------------------------------------------------------------- C SUMMARY OF FOUND TRACKS C----------------------------------------------------------------------- print 4011,n_mctrak_plausbl 4011 format(' tpc_hitlist:', 2 ' number of plausible tracks in list=',I4) if(n_mctrak_plausbl.gt.0)then do 320 i_mctrak_plausbl=1,n_mctrak_plausbl print 4012,i_mctrak_plausbl, 1 CU_mctrak_plausbl(i_mctrak_plausbl), 2 FI_mctrak_plausbl(i_mctrak_plausbl), 3 DZ_mctrak_plausbl(i_mctrak_plausbl) 4012 format(' ',I4, 2 ' C=',F10.5,' F=',F10.5,' DZ=',F10.5) 320 continue endif C end of condition on L_EVENT_MATCH and number of hits ENDIF C end of compliation flag CLEO_DOITTPC #endif RETURN END