#!/usr/bin/env python3 # ''' .. module:: HDF_Eventpar :synopsis: Pipeline script for X-ray event processing. .. moduleauthor:: Cor de Vries Suited for multi-processor execution using 'mpiexec'. ''' import numpy def mexit(rank,msg): if rank == 0: sys.exit(msg) else: sys.exit() def replicate(n,data): pdata=[] for i in numpy.arange(n): pdata.append(data) return pdata # ============================================================= if __name__ == "__main__": import sys,os import argparse parser = argparse.ArgumentParser(\ description="Process event HDF5 records for optimal filtering and spectral fitting",\ epilog="Execution can be distributed over multiple CPU's, using 'mpiexec'") parser.add_argument('-ch','--channel',type=int,required=False,default=None,\ help='Channel number to use. Default is first available channel.') parser.add_argument('-f','--frequency',type=str,nargs='+',required=False,default=[None],\ help="Frequency numbers (pixels) to use. Default is first available frequency. 'all' for all frequencies ") parser.add_argument('-n','--nrecs',type=str,required=False,default='all',\ help='Number of records to display. When single is not selected, average of records is shown. Default is all records') parser.add_argument('--delist',type=str,required=False,default=None,\ help='file specifying the event records numbers to ignore; list of "from-to" to ignore') parser.add_argument('-sum','--summary',action='store_true',required=False,\ help='Print summary of results to file') parser.add_argument('--pdf',action='store_true',required=False,\ help='Only pdf output') parser.add_argument('--debug',action='store_true',required=False,\ help='show some debug info + plots on screen') parser.add_argument('--threshold','-t',type=int,default=2000,required=False,\ help='Amplitude threshold for event detection (default 2000)') parser.add_argument('--bsec',type=float,required=False,default=0.05,\ help='Section (start+end) of record to use for background (default 0.05)') parser.add_argument('--bpo',action='store_true',required=False,\ help='select only pretrigger section for baseline value') parser.add_argument('--nopulseselect',action='store_true',required=False,\ help='discard event selection based on pulse intensity') parser.add_argument('--absolute',action='store_true',required=False,\ help='Use sqrt(I^2+Q^2) for pulse record') parser.add_argument('--noise',action='store_true',required=False,\ help='Let the script ask for seperate noise file') parser.add_argument('--brange',nargs=2,type=int,default=[None,None],\ help='Range of baseline values to select. Default: no selection') parser.add_argument('-nf','--noisefile',type=str,required=False,default=None,\ help='Define seperate noise file(name)') parser.add_argument('-np','--nopulseposition',action='store_true',required=False,\ help='do not select pulses based on pulse position') parser.add_argument('-nc','--nocorrelation',action='store_true',required=False,\ help='do not correct for opt.filter vs. baseline correlation') parser.add_argument('-fc','--freqcutoff',type=int,required=False,default=None,\ help='Maximum frequency cutoff for optimal filtering') parser.add_argument('-pl','--prlen',type=int,required=False,default=None,\ help='Length (power of 2) of record to process for optimal filtering (default entire record)') parser.add_argument('-rot','--rotate',action='store_true',required=False,\ help='Rotate events for minimum Q-signal (based on 100 events phase average)') parser.add_argument('-ph','--phase',type=float,required=False,default=None,\ help='Rotate events with given fixed phase (degrees)') parser.add_argument('--tdfilter',action='store_true',required=False,\ help='Perform optimal filtering in time domain') parser.add_argument('-ect','--ectfilter',action='store_true',required=False,\ help='Perform filtering to exclude electric crosstalk positive pulses') parser.add_argument('-rt','--risetime',action='store_true',required=False,\ help='compute and show histograms of rise time for pulses (only for single frequency)') parser.add_argument('--rterange',type=float,nargs=2,default=[None,None],\ help='range in energy for rise time correlation plot') parser.add_argument('--afitrange',type=float,nargs=2,default=[5.860,5.920],\ help='Fitrange in keV for alpha lines (default 5.860 5.920 )') parser.add_argument('--avout',action='store_true',required=False,\ help='output ASCII text files with average pulses') parser.add_argument('-ss','--sumspectrum',action='store_true',required=False,\ help='show spectrum, summed over all frequencies (pixels)') parser.add_argument('-el','--eventlist',action='store_true',required=False,\ help='write eventlist with optimal filter and baseline values') parser.add_argument('-rel','--raweventlist',action='store_true',required=False,\ help='write eventlist of all pulses with raw pulse parameters (no optimal fit and energy)') parser.add_argument('--eventtime',action='store_true',required=False,\ help='write event times in eventlist file') parser.add_argument('--woptf',action='store_true',required=False,\ help='write optimal filter data to file') parser.add_argument('--useoptf',type=str,required=False,default=None,\ help='use given external optimal filter file template') parser.add_argument('--etprev',type=float,nargs=2,default=[None,None],required=False,\ help='only select events with delta time within given interval from previous event') parser.add_argument('--etnext',type=float,nargs=2,default=[None,None],required=False,\ help='only select events with delta time within given interval to next event') parser.add_argument('--alltim',action='store_true',required=False,\ help='for deltatimes for "etprev", use all events from all pixels') parser.add_argument('dir',help='Directory or file of x-ray data HDF5 file(s)') from mpi4py import MPI comm=MPI.COMM_WORLD # start multiprocess administration rank=comm.rank # identify this process nproc=comm.size # total number of processes if rank == 0: # use parser only in process 0, and scatter results args=parser.parse_args() if args.debug: args.pdf=False pargs=replicate(nproc,args) print("Number processes used: ",nproc) else: if ( len(sys.argv) > 1 ) and ( sys.argv[1] == '-h'): sys.exit() pargs=None args=comm.scatter(pargs,root=0) freqcutoff=args.freqcutoff import matplotlib if args.pdf: matplotlib.use('Agg') else: matplotlib.use('GTK3Agg') import matplotlib.pyplot as plt from matplotlib.backends.backend_pdf import PdfPages from tesfdmtools.utils.widgets.filechooser import getfile from tesfdmtools.hdf.HMUX import HMUX from tesfdmtools.methods.fitcor import fitcor from tesfdmtools.methods.raweventlist import raweventlist_out from tesfdmtools.methods.positivepulse import positivepulse from tesfdmtools.utils.fiterrplot import fiterrplot if args.tdfilter: from tesfdmtools.methods.HDF_Eventutils import \ eventpars,plotepar,selecthist,selectevs,noisefile,noisespec,noisplot,getphase,\ pulseshape,getrisetime,limdist,showcorr,timplt,crosscplt,makespectrum,\ fitspectrum,makeplot,xoutl,plotptimes,eventlist_out,subsel,plotepar,selectdtim if args.woptf or ( args.useoptf is not None): from tesfdmtools.methods.TDMoptfilter import optfilter else: from tesfdmtools.methods.TDoptfilter import optfilter else: if args.woptf or ( args.useoptf is not None): from tesfdmtools.methods.HDF_Eventutils import \ eventpars,plotepar,selecthist,selectevs,noisefile,noisespec,noisplot,getphase,\ pulseshape,getrisetime,limdist,showcorr,timplt,crosscplt,makespectrum,\ fitspectrum,makeplot,xoutl,plotptimes,eventlist_out,subsel,plotepar,selectdtim from tesfdmtools.methods.TMoptfilter import optfilter else: from tesfdmtools.methods.HDF_Eventutils import \ eventpars,plotepar,selecthist,selectevs,noisefile,noisespec,noisplot,getphase,\ pulseshape,getrisetime,optfilter,limdist,showcorr,timplt,crosscplt,makespectrum,\ fitspectrum,makeplot,xoutl,plotptimes,eventlist_out,subsel,plotepar,selectdtim from tesfdmtools.methods.holzer import holzer,holzerfit from tesfdmtools.methods.filterect import filterect if os.path.isfile(args.dir): cfilename=args.dir else: if rank == 0: # only select filename in process 0 if not os.path.exists(args.dir): print('Error, directory %s does not exist.' % args.dir) cfilename=getfile(pattern='*.h5',path=args.dir) pcfilename=replicate(nproc,cfilename) else: pcfilename=None cfilename=comm.scatter(pcfilename,root=0) # scatter selected filename to all processes if args.prlen is not None : # optimal filter procesing length proclen=2**int(numpy.log(args.prlen)/numpy.log(2)) # must be power of 2 else: proclen=None hdf=HMUX(filename=cfilename) fnam=os.path.basename(cfilename).split('.')[0] channel=args.channel if channel is None: channel=int(hdf.channels[0]) conf_tab=hdf.channel[channel].freq_conf cpixels=conf_tab['pixel_index'] fpixels=hdf.channel[channel].freqs apixels=numpy.intersect1d(cpixels,fpixels) # take only frequencies which occur both in the table as in the channel if rank == 0: print("Available pixels(frequencies): ",apixels) if apixels.size == 0: mexit(rank,"pixel configuration table and available frequencies do not match") if args.frequency[0] is None: frequencies=[apixels[0]] else: fmt='Error, illegal frequency specification: "'+len(args.frequency)*'%s '+'"' if args.frequency[0] == 'all': if len(args.frequency) > 1: mexit(rank,fmt % tuple(args.frequency)) else: frequencies=numpy.array(sorted(apixels),dtype=int) else: frequencies=[] for fff in args.frequency: if fff.isdigit(): frr=int(fff) findx,=numpy.where(apixels == frr) if len(findx) == 0: mexit(rank,"Selected frequency (pixel) %d does not exist in file" % frr) frequencies.append(frr) else: mexit(rank,fmt % tuple(args.frequency)) frequencies=numpy.array(sorted(frequencies),dtype=int) nrow=4 ncol=1 fontsize=10 if len(frequencies) > 4: ncol=2 fontsize=7 else: nrow=len(frequencies) pid=os.getpid() # make unique id for this process host=os.uname()[1] unq="%s_%d_" % (host,pid) matplotlib.rcParams.update({'font.size': fontsize}) show='X' if args.pdf: show='pdf' rotate=args.rotate rotp=0.0 if args.phase is not None: rotate=args.phase/180.0*numpy.pi rotp=rotate ptitle='%s' % fnam pdfplt={} fitplt={} crsplt={} nosplt={} titplt={} nfile=None # if len(frequencies) == 1: # risetime plots only for single frequency # risetime=args.risetime # else: # risetime=False risetime=args.risetime if rank == 0: # identify info for the separate processes if args.summary: sfile=open('Eventpar_'+fnam+'.txt','w') sfile.write('file: %s\n' % fnam) sfile.write('pars:') pttt='' for pp in sys.argv[1:-1]: pttt=pttt+' '+pp pttt=pttt+'\n' sfile.write(pttt) sfile.write('\n i px freq(khz) V_bias gbwp fwhm(err)(eV) counts alin.fct rise fall\n\n') tdata={} # list of iplt's and frequencies to spread over processes for iproc in numpy.arange(nproc): # initialize tdata[iproc]=[] iproc=0 for iplt,freq in enumerate(sorted(frequencies)): # spread frequencies to be processed over different processes if iproc in tdata: tdata[iproc].append([iplt,freq]) else: tdata[iproc]=[[iplt,freq]] iproc=iproc+1 if iproc >= nproc: iproc=0 pdata=[] for pp in tdata: pdata.append(tdata[pp]) else: pdata=None frequencies=comm.scatter(pdata,root=0) # now scatter info to all processes if args.alltim: # look at time of all event from all pixels from tesfdmtools.methods.aldtimes import aldtimes,selectald aldtims=aldtimes(comm,rank,hdf,channel=channel,frequencies=frequencies,nrecs=args.nrecs,\ threshold=args.threshold,debug=args.debug) sfiletxt={} # process events into spectra for iplt,freq in frequencies: print("process: ",unq," (iplt,freq): ",iplt,freq) pdfnam='freq%3.3d_Events_%s.pdf' % (iplt,fnam) # unique plot file for each frequency pdf = PdfPages(pdfnam) pdfplt[iplt]=pdfnam pindx,=numpy.where(conf_tab['pixel_index'] == freq) if ( type(rotate) is bool ) and ( rotp == 0.0 ): # any phase rotation to be computed? if rotate: indx=numpy.zeros(100) ii=0 for j,record in enumerate(hdf.channel[channel].freq[freq]): # get index of first 100 valid records if numpy.ptp(record) != 0: indx[ii]=j ii=ii+1 if ii == 100: break rotp=getphase(hdf,channel,freq,indx,debug=args.debug) print("Computed phase rotation: ",rotp) fffx,=numpy.where( hdf.channel[channel].freq_conf['pixel_index'] == freq ) freqindx=fffx[0] khz=hdf.channel[channel].freq_conf['freq'][freqindx] flip=False if khz == 0.0: print("Data for channel %d, freq %d has zero frequency. Data are DC data" % (channel,freq)) flip=positivepulse(hdf.channel[channel].freq[freq],args.threshold) # flip data in record for DC # (frequency=0.0) positive pulse data channel,freq,blines1,blines2,positions,peaks,surf,ptp,evttim,rclen=\ eventpars(hdf,channel=channel,freq=freq,nrecs=args.nrecs,rotation=rotp,\ threshold=args.threshold,bsec=args.bsec, bpt=args.bpo,delist=args.delist,flip=flip) # get basic pulse parameters if args.debug: print ("freq: %d nevttime: %d\n" % (freq,evttim.size) ) np=20 pfmt="evttim: "+np*" %6.2f"+"\n" print (pfmt % tuple(evttim[0:np])) # bline=(blines1+blines2)/2.0 bline=blines1+positions/float(rclen)*(blines2-blines1) raweventlist_out(args.raweventlist,cfilename,freq,bline,positions,peaks,surf,evttim,channel=channel) indx=selectevs(hdf,channel,freq,bline,positions,peaks,surf,rclen,\ show=show,debug=args.debug,nopos=args.nopulseposition,\ pdf=pdf,nopulse=args.nopulseselect,bsec=args.bsec) # select proper events if args.alltim and ( args.etprev[0] is not None ): # select events based on prev time to any event indx=selectald(indx,aldtims[freq],etprev=args.etprev,debug=args.debug) else: if ( args.etprev[0] is not None ) or \ ( args.etnext[0] is not None ): # select events on time to prev/next event indx=selectdtim(indx,evttim,etprev=args.etprev,etnext=args.etnext) if ( args.brange[0] is not None ) and \ ( args.brange[1] is not None ): # any additional baseline selection? indx=subsel(indx,bline,args.brange) if args.ectfilter: ectpars=filterect(hdf,channel,freq,indx,blines1,blines2) # filter electric crosstalk positive pulses else: ectpars=None if show is not None: plotepar(hdf,channel,freq,bline[indx],positions[indx],peaks[indx],surf[indx],\ show=show,pdf=pdf,indx=indx,ectpar=ectpars) # show parameters of selected events if args.ectfilter: indx=indx[ectpars[0]] # use electric crosstalk filter if applicable if indx.size < 20: print("Insufficient events for frequency: ",freq) fitplt[iplt]=None titplt[iplt]="ch=%d, px=%d" % (channel,freq) # plot titles crsplt[iplt]=None nosplt[iplt]=None # os.unlink(pdfplt[iplt]) pdf.close() continue if args.noise: # seperate noise file ? nfile=True # ask for file else: nfile=args.noisefile # noise file can be given nhdf,nptp=noisefile(hdf,channel,freq,nfile=nfile) if nhdf is None: noisespc,fax=noisespec(hdf,channel,freq,ptp,debug=args.debug,\ absolute=args.absolute,prlen=proclen,flip=flip) # compute noise spectrum from xray file else: noisespc,fax=noisespec(nhdf,channel,freq,nptp,debug=args.debug,\ absolute=args.absolute,prlen=proclen,flip=flip) # compute noise spectrum from noise file avpulsepos=numpy.mean(positions[indx]) # average pulseposition ifit,rtimes,ftimes,avpulse,avbline=optfilter(hdf,channel,freq,indx,blines1,blines2,noisespc,\ freqcutoff=freqcutoff,debug=args.debug,rotate=rotate,\ risetime=risetime,bsec=args.bsec,absolute=args.absolute,\ prlen=proclen,apulsepos=avpulsepos,\ wrtfilter=args.woptf,usefilter=args.useoptf,flip=flip) # compute optimal filter fits avspc=numpy.absolute(numpy.fft.fft(avpulse)[0:(avpulse.size//2)]) nosplt[iplt]={'nspectrum':noisespc,'spectrum':avspc,'fax':fax} if args.nocorrelation: nfit=ifit center=None ytilt=None else: nfit,center,ytilt=fitcor(bline[indx],ifit,debug=args.debug) # compute optimal-filter-fit versus baseline correlation if args.debug and ( not args.nocorrelation): showcorr(bline[indx],ifit,xlabel='baseline',ylabel='e-fit',\ title=('cross_corr, px=%d' % freq),center=center,ytilt=ytilt) evtttt=None if args.eventtime: evtttt=evttim[indx] eventlist_out(args.eventlist,freq,indx,nfit,bline[indx],peaks[indx],cfilename,\ channel=channel,rtimes=rtimes,ftimes=ftimes,evttim=evtttt,surf=surf[indx]) # write eventlist spectrum,bins=makespectrum(nfit,debug=args.debug) # make spectrum crsplt[iplt]={'ifit':ifit,'bline':bline[indx],'center':center,'ytilt':ytilt,\ 'frange':[bins[0],bins[-1]],'index':indx} # store correlation parameters fitplt[iplt]={} fwhm,fpars=fitspectrum(fitplt,iplt,hdf,channel,freq,spectrum,bins,\ debug=args.debug,afitrange=args.afitrange) # fit spectral lines titplt[iplt]="ch=%d, px=%d %7.2f (kHz)" % (channel,freq,fpars['khz']) # plot titles fwhmtxt='%8.2f(%4.2f)' % (fwhm,fpars['fwhmerr']) # fwhm + error (eV) if risetime: pttitle="%s ch=%d, px=%d" % (ptitle,channel,freq) plotptimes(rtimes,ftimes,nfit,avpulse,pttitle,pdf,ps=args.pdf,outtxt=args.avout,\ erange=args.rterange,avbline=avbline) # risetime plot xix=xoutl(rtimes) avrt=numpy.mean(rtimes[xix])*1E6 xix=xoutl(ftimes) avft=numpy.mean(ftimes[xix])*1E6 rtxt='%6.1f %6.1f' % (avrt,avft) else: rtxt='' # store information of this process if args.summary: sfiletxt[freq]='%3d %3d %11.4f %9.2f %9.4f %14.14s %8d %8.4f %s\n' % \ (iplt,freq,fpars['khz'],conf_tab['ampl'][pindx],conf_tab['gbwp'][pindx],\ fwhmtxt,fpars['counts'],fpars['lratio'],rtxt) pdf.close() # store information of this process sdata={'sfiletxt':sfiletxt,'pdfnam':pdfplt,'fitplt':fitplt,'crsplt':crsplt,\ 'nosplt':nosplt,'titplt':titplt,'frequencies':frequencies} # gather information from all processes pdata=comm.gather(sdata,root=0) if rank != 0: # stop all processes except main process print("Stop process %s" % unq) sys.exit() # for the main process start assembling of the data print("finishing up in process %s...................." % unq) pdfs={} sfiletxt={} fitplt={} crsplt={} nosplt={} titplt={} fff=[] for data in pdata: for iplt in data['pdfnam']: pdfs[iplt]=data['pdfnam'][iplt] for ff in data['sfiletxt']: sfiletxt[ff]=data['sfiletxt'][ff] for iplt in data['fitplt']: fitplt[iplt]=data['fitplt'][iplt] for iplt in data['crsplt']: crsplt[iplt]=data['crsplt'][iplt] for iplt in data['nosplt']: nosplt[iplt]=data['nosplt'][iplt] for iplt in data['titplt']: titplt[iplt]=data['titplt'][iplt] for iplt,frq in data['frequencies']: fff.append(frq) frequencies=numpy.array(sorted(fff),dtype=int) pdfmain='main_Events_%s.pdf' % fnam pdf = PdfPages(pdfmain) # line fit plots for jplt in numpy.arange(len(frequencies)): iplt=jplt % 8 if iplt == 0: plots={} plt.suptitle(ptitle,size=12) plots[iplt] = plt.subplot2grid((nrow,ncol), ((iplt//ncol),(iplt % ncol)) ) fp=fitplt[jplt] if fp is not None: fiterrplot(fp['xax'],fp['spectrum'],fp['yerr'],fp['holy'],\ xtitle="Energy (keV)",ytitle="Counts",ptitle=fp['ptitle'],ptxt=fp['ptxt'],pltax=plots[iplt]) else: plots[iplt].set_title(titplt[jplt]) if ( iplt == 7 ) or ( jplt == (len(frequencies)-1) ): # psfile="Events_lif_%s.ps" % fnam makeplot(pdf,ps=args.pdf) # summed spectrum ( total of all spectra for all pixels ) if args.sumspectrum and ( len(frequencies) > 1 ): for j,iplt in enumerate(fitplt): fp=fitplt[iplt] if fp is not None: if j == 0: xax=fp['xax'] tspec=fp['spectrum'] tspec[0]=0 tspec[-1]=0 else: ss=fp['spectrum'] ss[0]=0 ss[-1]=0 tspec=tspec+numpy.interp(xax,fp['xax'],ss) wsspect=True if wsspect: wsnam='Summed_spectrum_%s.txt' % fnam wsfile=open(wsnam,'w') wsfile.write('# Summed spectrum: %s\n' % fnam) wsfile.write('# Energy Counts\n') for j,xx in enumerate(xax): wsfile.write('%12.6f %12d\n' % (xx,tspec.astype(int)[j]) ) wsfile.close() ppp=holzerfit(xax,tspec.astype(int)) # fit total spectrum holy=holzer(ppp,xax) # fitted function yrr=numpy.sqrt(tspec) # 1-sigma errors yerr=numpy.where(yrr <= 0.0,1.0,yrr) # for zero counts insert one-count error fwhmtxt="%5.5s %12.4g +/- %12.4g (eV)" % ("FWHM: ",(ppp[3]*1000.0),(ppp[9]*1000.0)) # FWHM + 1-sigma error sttt="Average resolution for all pixels %5.5s %12.4g +/- %12.4g (eV)" % \ ("FWHM: ",(ppp[3]*1000.0),(ppp[9]*1000.0)) sfiletxt['total']='\n\n'+' '.join(sttt.split()) sfiletxt['totct']='\nTotal counts: %d' % (numpy.sum(tspec)) ptxt=' '.join(fwhmtxt.split())+'\n' ptxt=ptxt+( "%d counts\n" % (numpy.sum(tspec)) ) ptxt=ptxt+( "C-stat: %5.1f\n" % ppp[4] ) ptxt=ptxt+( "degs of frdm: %d\n" % ppp[5] ) plt.suptitle(ptitle,size=12) splot=plt.subplot2grid((2,1), (0,0) ) fiterrplot(xax,tspec,yerr,holy,\ xtitle="Energy (keV)",ytitle="Counts",ptitle="Summed spectrum",ptxt=ptxt,pltax=splot) makeplot(pdf,ps=args.pdf) # (ascii) summary file if args.summary: for ff in frequencies: if ff in sfiletxt: sfile.write(sfiletxt[ff]) else: print("No summary info found for frequency: ",ff) if 'total' in sfiletxt: sfile.write(sfiletxt['total']) sfile.write(sfiletxt['totct']) sfile.write('\n') sfile.close() # time evolution plots for jplt,freq in enumerate(frequencies): iplt=jplt % 8 if iplt == 0: plots={} plt.suptitle(ptitle,size=12) plots[iplt] = plt.subplot2grid((nrow,ncol), ((iplt//ncol),(iplt % ncol)) ) fp=crsplt[jplt] if fp is not None: timplt(plots[iplt],hdf,channel,freq,fp['index'],fp['ifit'],frange=fp['frange'],\ xlabel='record number',ylabel='e-fit') if ( iplt == 7 ) or ( jplt == (len(frequencies)-1) ): # psfile="Events_crs_%s.ps" % fnam makeplot(pdf,ps=args.pdf) # crosscorr plots if not args.nocorrelation: for jplt,freq in enumerate(frequencies): iplt=jplt % 8 if iplt == 0: plots={} plt.suptitle(ptitle,size=12) plots[iplt] = plt.subplot2grid((nrow,ncol), ((iplt//ncol),(iplt % ncol)) ) fp=crsplt[jplt] if fp is not None: crosscplt(plots[iplt],hdf,channel,freq,fp['bline'],fp['ifit'],\ xlabel='baseline',ylabel='e-fit',center=fp['center'],ytilt=fp['ytilt']) if ( iplt == 7 ) or ( jplt == (len(frequencies)-1) ): # psfile="Events_crs_%s.ps" % fnam makeplot(pdf,ps=args.pdf) # noise spectrum plots for jplt,freq in enumerate(frequencies): iplt=jplt % 8 if iplt == 0: plots={} plt.suptitle(ptitle,size=12) plots[iplt] = plt.subplot2grid((nrow,ncol), ((iplt//ncol),(iplt % ncol)) ) fp=nosplt[jplt] if fp is not None: title='ch=%d px=%d %7.2f (kHz)' % (channel,freq,conf_tab['freq'][iplt]) noisplot(plots[iplt],fp['spectrum'],fp['nspectrum'],fp['fax'],title) if ( iplt == 7 ) or ( jplt == (len(frequencies)-1) ): # psfile="Events_nos_%s.ps" % fnam makeplot(pdf,ps=args.pdf) pdf.close() # now combine pdf's and delete duplicates print("combine pdf's................") cmd='pdfunite %s' % pdfmain for iplt in numpy.arange(len(frequencies)): if os.path.exists(pdfs[iplt]): cmd=cmd+' '+pdfs[iplt] pdfnam='Events_%s.pdf' % fnam cmd=cmd+' '+pdfnam # print(cmd) xx=os.system(cmd) if xx == 0: os.unlink(pdfmain) for pnn in pdfs: if os.path.exists(pdfs[pnn]): os.unlink(pdfs[pnn]) else: sys.exit("Error executing 'pdfunite' to combine pdf's ") print("Finished!")