c Searches though BATSE time data, looks for disconnected echo c of same hardness. c Version 2.0: version running through many GRBs. c Version 3.0: searches for peak in summed 4 channel data. c Version 3.1: variance of total counts now used in NR chi**2. c Version 3.2: background taken as last iwide bin. program millisearch implicit real*8 (a-z) integer*4 i, it, iw, itmax, ibin, itrig, npts, nlasc, preb integer*4 ie, iwide, iwidemax, iwmax, iwtop, idatagap integer*4 imax_burst, iburst_start, iburst_end integer*4 imax_echo, iecho_start, iecho_end integer*4 imb character atrig*5, prefix*32, fname*37, acomment*30 parameter (itmax=100000) real*8 c1(itmax), c2(itmax), c3(itmax), c4(itmax), cT(itmax) real*8 w1(itmax), w2(itmax), w3(itmax), w4(itmax), wT(itmax) real*8 burst(4), burst_back(4), echo(4), echo_back(4) real*8 Rburst(4), Rvariance(4), Secho(4), Svariance(4) c Open files. open (1, file='/phyfac23/batse/cat64ms/triglist.in', 2 status='old') open (2, file='millisearch.out', status='unknown') write (2,*) ' itrig, fmax, tmin, tmax, chi2, acomment' c do 800 igrb=1, 10000 read (1,2,end=999) atrig 2 format (54x, a5) c Zero out needed arrays. do 4 i=1, itmax c1(i)=0.0d0 c2(i)=0.0d0 c3(i)=0.0d0 c4(i)=0.0d0 cT(i)=0.0d0 w1(i)=0.0d0 w2(i)=0.0d0 w3(i)=0.0d0 w4(i)=0.0d0 wT(i)=0.0d0 4 continue c itrig=int(atrig) prefix='/phyfac23/batse/cat64ms/cat64ms.' fname=prefix//atrig write (*,*) fname c open (3, file=fname, status='old') c Read in header information. read (3,*) read (3,*) itrig, npts_in, nlasc, preb c Read in flux data. c write (*,*) ' Reading in data.' do 10 i=1, npts_in read (3,*,end=15) c1(i), c2(i), c3(i), c4(i) cT(i)=c1(i) + c2(i) + c3(i) + c4(i) 10 continue 15 npts=i-1 c Check for data in file. acomment=' No data in file.' cprod=c1(1)*c2(1)*c3(1)*c4(1) if (cprod.le.1.0d-6) goto 780 c Check to see if any data gaps. c write (*,*) ' Checking for data gaps.' idatagap=0 acomment=' Data gap detected.' do 20 i=1, npts cprod=c1(i)*c2(i)*c3(i)*c4(i) if (cprod.le.1.0d-6) idatagap=1 20 continue c End data at beginning of data gap. if (idatagap.eq.1) npts=i-1 c Check slope of background. c See if any bins drop below backsig c of the initial bin. c write (*,*) ' Checking for sloping, curving background.' acomment=' Background too hostile.' backsig=4.5d0 s1=c1(1)-backsig*sqrt(c1(1)) s2=c2(1)-backsig*sqrt(c2(1)) s3=c3(1)-backsig*sqrt(c3(1)) s4=c4(1)-backsig*sqrt(c4(1)) do 30 i=1, npts if (c1(i).lt.s1) goto 780 if (c2(i).lt.s2) goto 780 if (c3(i).lt.s3) goto 780 if (c4(i).lt.s4) goto 780 30 continue c c************************************************************* c c This GRB has no disqualifying data gaps or background. c Now find the highest S/N series of bins. c Call this the "burst." c c************************************************************* c c Find highest flux bin. c write (*,*) ' Finding highest flux bin.' cmax=0.0d0 do 100 i=1, npts-1 if (cT(i).le.cmax) goto 100 cmax=cT(i) imax_burst=i 100 continue c Find S/N for this/next bin. c write (*,*) ' Finding best S/N binning around high.' imb=imax_burst cmain=cT(imb) cback=cT(imb+1) iwide=1 iburst_start=imax_burst cmainmax=cmain cbackmax=cback signal=cmain-cback acomment=' Negative signal.' if (signal.le.-1.0d-6) goto 780 noise=sqrt(cback) SdN=signal/noise SdNmax=SdN c Add bins RIGHT to find best S/N. do 200 it=1, (npts-imax_burst-1)/2 cmain=cmain+cT(imax_burst+it) cback=cback-cT(imax_burst+it)+cT(imax_burst+it+it) 2 +cT(imax_burst+it+it+1) c Is moving back withing 2sig of final back? c ncompare=abs(cback/float(it+1)-cT(npts)) c if (ncompare.gt.2.0d0*sqrt(cback)) goto 200 c Compute signal. signal=cmain-cback if (signal.le.1.0d-6) goto 200 noise=sqrt(cback) SdN=signal/noise c Is S/N better than ever? if (SdN.lt.SdNmax) goto 200 SdNmax=SdN iwide=it+1 iburst_start=imax_burst cmainmax=cmain cbackmax=cback 200 continue c Now add bins LEFT if increases S/N. cmain=cmainmax cback=cbackmax iwidemax=iwide c do 300 it=1, npts if (imax_burst-it.lt.1) goto 300 if (imax_burst+iwide+iwide+it-1.gt.npts) goto 300 cmain=cmain + cT(imax_burst-it) cback=cback + cT(imax_burst+iwide+iwide+it-1) c ncompare=abs(cback/float(iwide+it)-cT(npts)) c if (ncompare.gt.2.0d0*sqrt(cback)) goto 300 c Compute signal. signal=cmain-cback if (signal.le.1.0d-6) goto 300 noise=sqrt(cback) SdN=signal/noise c Is S/N better than ever? if (SdN.lt.SdNmax) goto 300 SdNmax=SdN iwidemax=iwide+it iburst_start=imax_burst-it cmainmax=cmain cbackmax=cback 300 continue c************************************************************* c c Rebin around S/N peak with new wider time bins. c c************************************************************* c write (*,*) ' Rebinning with new wider time bins.' iwide=iwidemax iburst_end=iburst_start + iwide - 1 iwmax=npts/iwide do 400 iw=1, iwmax if (iburst_start-1+iw*iwide.gt.npts) goto 450 do 350 it=1, iwide ibin=iburst_start+(iw-1)*iwide+(it-1) acomment=' iburst_start overwrite.' if (ibin.gt.npts) goto 780 w1(iw)=w1(iw) + c1(ibin) w2(iw)=w2(iw) + c2(ibin) w3(iw)=w3(iw) + c3(ibin) w4(iw)=w4(iw) + c4(ibin) wT(iw)=wT(iw) + cT(ibin) 350 continue 400 continue 450 iwtop=iw-1 c Reassign burst variables. burst(1)=w1(1) burst(2)=w2(1) burst(3)=w3(1) burst(4)=w4(1) burstT=burst(1)+burst(2)+burst(3)+burst(4) burst_back(1)=w1(iwtop) burst_back(2)=w2(iwtop) burst_back(3)=w3(iwtop) burst_back(4)=w4(iwtop) burst_backT=w1(iwtop)+w2(iwtop)+w3(iwtop)+w4(iwtop) c Is burst bright enough for analysis? acomment=' GRB not 5.5sig over backgd.' if (burstT-burst_backT.le.5.5d0*sqrt(burst_backT)) 2 goto 780 c************************************************************* c c Find second highest peak in rebinned wide data. c Call this the "echo." c c************************************************************* c write (*,*) ' Finding second highest time bin.' cTecho=cT(1) imax_echo=1 c do 500 it=imax_burst+iwide+iwide, npts if (cT(it).le.cTecho) goto 500 cTecho=cT(it) imax_echo=it 500 continue c acomment=' cT(1) > any echo bin.' if (imax_echo.eq.1) goto 780 iecho_start=imax_echo-(imax_burst-iburst_start) iecho_end=iecho_start + iwide - 1 acomment=' Echo overflows data start.' if (iecho_start-iwide+1.lt.1) goto 780 acomment=' Echo overflows data end.' if (iecho_end.gt.npts) goto 780 c echo(1)=0.0d0 echo(2)=0.0d0 echo(3)=0.0d0 echo(4)=0.0d0 echo_back(1)=0.0d0 echo_back(2)=0.0d0 echo_back(3)=0.0d0 echo_back(4)=0.0d0 c do 550 it=iecho_start, iecho_end echo(1)=echo(1) + c1(it) echo(2)=echo(2) + c2(it) echo(3)=echo(3) + c3(it) echo(4)=echo(4) + c4(it) echo_back(1)=echo_back(1) + c1(it-iwide) echo_back(2)=echo_back(2) + c2(it-iwide) echo_back(3)=echo_back(3) + c3(it-iwide) echo_back(4)=echo_back(4) + c4(it-iwide) 550 continue c Record tmin and tmax times. tmin=2.0d0*dble(iwide)*0.064d0 tmax=dble(npts-iwide+1-iburst_start)*0.064d0 acomment=' tmax .lt. tmin.' if (tmax.lt.tmin) goto 780 c Is echo some sig over background? echoT=echo(1) + echo(2) + echo(3) + echo(4) echo_backT=echo_back(1) + echo_back(2) 2 + echo_back(3) + echo_back(4) ereject=8.0 c if (echoT-echo_backT.lt.5.5d0*sqrt(echo_backT)) then if (echoT-echo_backT.lt.ereject*sqrt(echo_backT)) then c Too weak: record fmin echo strength. fmax=ereject*sqrt(echo_backT)/(burstT-burst_backT) chi2=1000.0d0 goto 650 else fmax=(echoT-echo_backT)/(burstT-burst_backT) endif c************************************************************* c c See if burst and echo have the same hardness c to within statistical errors. Use 2 bin chi**2. c See Numerical Recipes, 2nd edition, top of page 617. c c************************************************************* do 575 ie=1, 4 Rburst(ie)=burst(ie)-echo_back(ie) if (Rburst(ie).lt.0.0d0) Rburst(ie)=0.0d0 Rvariance(ie)=burst(ie) Secho(ie)=echo(ie)-echo_back(ie) if (Secho(ie).lt.0.0d0) Secho(ie)=0.0d0 Svariance(ie)=echo(ie) 575 continue c Rtot=Rburst(1)+Rburst(2)+Rburst(3)+Rburst(4) Stot=Secho(1)+Secho(2)+Secho(3)+Secho(4) chi2=0.0d0 c do 600 ie=1, 4 top=sqrt(Stot/Rtot)*Rburst(ie)-sqrt(Rtot/Stot)*Secho(ie) bot=Rvariance(ie) + Svariance(ie) chi2=chi2 + top**2/bot 600 continue c 650 acomment=' OK' c write (*,*) ' npts, imax_burst, imax_echo = ', c 2 npts, imax_burst, imax_echo c write (*,*) ' c3(imax_burst), c3(imax_echo) = ', c 2 c3(imax_burst), c3(imax_echo) c write (*,*) ' iburst_start, iburst_end = ', c 2 iburst_start, iburst_end c write (*,*) ' iecho_start, iecho_end = ', c 2 iecho_start, iecho_end c write (*,*) ' tmin, tmax = ', tmin/0.064, tmax/0.064 c write (*,*) ' fmax = ', fmax c write (*,*) ' chi**2 = ', chi2 c write (*,*) (w3(iw), iw=1, iwtop) c write (2,*) itrig, iburst_start, iburst_end, c 2 iecho_start, iecho_end c660 format (5(i8, 2x)) c 780 if (acomment.ne.' OK') then fmax=0.0d0 tmin=0.0d0 tmax=0.0d0 chi2=0.0d0 endif c write (2,785) itrig, fmax, tmin, tmax, chi2, acomment write (*,785) itrig, fmax, tmin, tmax, chi2, acomment 785 format (1x, i5, 2x, 3(f8.4, 2x), f8.1, 2x, a30) 800 continue c 999 stop end