#!/usr/bin/env python #./search_var.py M31Core PN500_10000_ 30 4 10 #./search_var.py M31Core PN500_10000_ 50 3 15 import sys import os import subprocess import glob import math import pyfits def selectEvents(eventfile,ra,dec,radius,outroot,iradius,oradius,boutroot): #write ds9 source region with open('%s.reg' %(outroot),'w+') as f: f.write('fk5\n') f.write('circle(%s,%s,%s") #text={%s}\n' %(ra,dec,radius,outroot)) pass #write ds9 background region with open('%s.reg' %(boutroot),'w+') as f: f.write('fk5\n') f.write('annulus(%s,%s,%s",%s") #text={%s}\n' %(ra,dec,iradius,oradius,boutroot)) pass #Run fselect filtering the event file with the regions we just created cmdLine = '''fselect infile=%s+1 outfile=%s.ds expr="regfilter('%s.reg')" clobber=yes copyall=no ''' % (eventfile,outroot,outroot) print cmdLine subprocess.call(cmdLine,shell=True) cmdLine = '''fselect infile=%s+1 outfile=%s.ds expr="regfilter('%s.reg')" clobber=yes copyall=no ''' % (eventfile,boutroot,boutroot) print cmdLine subprocess.call(cmdLine,shell=True) newdir = sys.argv[1] root = sys.argv[2] trans_threshold = float(sys.argv[3]) var_threshold = float(sys.argv[4]) rat_threshold = float(sys.argv[5]) filtered = sys.argv[6] ext_threshold = float(sys.argv[7]) print('Newdir = %s, root = %s, trans_threshold = %s, var_threshold = %s, rat_threshold = %s ext_threshold = %s' %(newdir,root,trans_threshold,var_threshold,rat_threshold,ext_threshold)) os.chdir(newdir) fileList = glob.glob('%s*.out' %(root)) ptransients = [] pvariableSources = [] transients = [] variableSources = [] tsingle = [] vsingle = [] for myslice in fileList: with open(myslice,'r') as f: lines = f.readlines() pass for line in lines: elementi = line.split() ra1 = elementi[0] dec1 = elementi[1] transient = False variable = True try: ra1 = float(ra1) if(ra1==0.0): #Error continue else: #Variable source ra = ra1 dec = float(dec1) transient = False pass except: if(ra1=='sorgente'): #Transient ra = float(elementi[6]) dec = float(elementi[7]) transient = True pass pass #Compute significance if(transient): #Control that the source is not extended ext1 = elementi[13] # if(float(ext1)!=0): if(float(ext1) > ext_threshold): #The source is extended. Cannot be a transient, skip it continue pass #TBD sig = float(elementi[14]) if(sig >= trans_threshold): pass if(transient): transFeatures = " ".join(elementi[6:]) ptransients.append('%s: %s' %(myslice,transFeatures)) tsingle.append(True) pass pass else: #Control if the source is extended, either in the whole observation #or in the slice ext1 = elementi[15] ext2 = elementi[16] # if(float(ext1)!=0 or float(ext2)!=0): if(float(ext1) > ext_threshold or float(ext2) > ext_threshold): #The source is extended, skip it continue pass #Variable source rate1 = float(elementi[7]) rate1Err = float(elementi[8]) rate2 = float(elementi[9]) rate2Err = float(elementi[10]) rat = rate1/rate2 ratErr = rat * math.sqrt(pow(rate1Err/rate1,2)+pow(rate2Err/rate2,2)) sig = (rate1-rate2)/math.sqrt(pow(rate1Err,2)+pow(rate2Err,2)) if(abs(sig) >= var_threshold and (rat >= rat_threshold or rat <= 1./rat_threshold)): pass if(variable): srcFeatures = " ".join(elementi[0:]) pvariableSources.append('%s: %s %s %s' %(myslice,srcFeatures,rat,ratErr)) vsingle.append(True) pass pass pass pass pass for cvar,var in enumerate(pvariableSources): timeint = var.split('_')[2].split('.out')[0] print('timeint=%s' %(timeint)) ndet = 1 elements = var.split() ra,dec,poserr = (float(elements[1]),float(elements[2]),float(elements[3])) print('vsingle=%s' %(vsingle[cvar])) if (vsingle[cvar]): for covar,ovar in enumerate(pvariableSources): if (covar > cvar and vsingle[covar]): vtimeint = ovar.split('_')[2].split('.out')[0] ele = ovar.split() vra = float(ele[1]) vdec = float(ele[2]) verr = float(ele[3]) # dist = 3600.*math.sqrt(pow(vra-ra,2)+pow(vdec-dec,2)) dist = 3600.*math.sqrt(((pow(vra-ra,2))*(pow(math.cos(vdec/(57.29577951308232)),2)))+(pow(vdec-dec,2))) diste = 2.+math.sqrt(pow(verr,2)+pow(poserr,2)) print('RA1 = %s, RA2 = %s, dist=%s, disterr=%s, ndet=%s' %(ra,vra,dist,diste,ndet)) if (dist/diste <= 3): timeint=timeint+','+vtimeint ndet = ndet+1 vsingle[covar] = False pass variableSources.append('%s %s %s' %(var,ndet,timeint)) pass pass #print('NOTSorted:%s' %(ptransients)) def my_order(element): elements = element.split(":")[-1] # print("Elements is %s" % elements) tokens = elements.split() # print("Tokens is %s" % tokens) # print("Tokens[-3] is %s" % tokens[-3]) return float(tokens[-3]) sptransients=sorted(ptransients,key=my_order)[::-1] #print('Sorted:%s' %(sptransients)) for ctrans,trans in enumerate(sptransients): timeint = trans.split('_')[2].split('.out')[0] print('timeint=%s' %(timeint)) ndet = 1 elements = trans.split() ra,dec,poserr = (float(elements[1]),float(elements[2]),float(elements[3])) print('tsingle=%s' %(tsingle[ctrans])) if (tsingle[ctrans]): for cotrans,otrans in enumerate(sptransients): if (cotrans > ctrans and tsingle[cotrans]): ttimeint = otrans.split('_')[2].split('.out')[0] ele = otrans.split() tra = float(ele[1]) tdec = float(ele[2]) terr = float(ele[3]) # dist=((sqrt((((ra1[i]-ra2[j])*(ra1[i]-ra2[j]))*cos(dec1[i]/(57.29577951308232))*cos(dec1[i]/(57.29577951308232)))+((dec1[i]-dec2[j])*(dec1[i]-dec2[j]))))); # dist = 3600.*math.sqrt(pow(tra-ra,2)+pow(tdec-dec,2)) dist = 3600.*math.sqrt(((pow(tra-ra,2))*(pow(math.cos(tdec/(57.29577951308232)),2)))+(pow(tdec-dec,2))) diste = 2.+math.sqrt(pow(terr,2)+pow(poserr,2)) print('RA1 = %s, RA2 = %s, dist=%s, disterr=%s, ndet=%s' %(ra,tra,dist,diste,ndet)) if (dist/diste <= 3): timeint=timeint+','+ttimeint ndet = ndet+1 tsingle[cotrans] = False pass transients.append('%s %s %s' %(trans,ndet,timeint)) pass pass f = open('%strans.txt' %(root),'w+') f.write('List of transient sources:\n') f.write('ra1,dec1,errpos1,rate1,errate1,cts1,errcts1,ext1,errext1,ml1,bkg1,idsrc1\n') for trans in transients: f.write(trans+'\n') pass f.close() print('Wrote %s transient sources' %(len(transients))) with open('%svar.txt' %(root),'w+') as f: f.write('List of variable sources:\n') f.write('ra1,dec1,errpos1,ra2,dec2,errpos2,dist,rate1,errate1,rate2,errate2,cts1,errcts1,cts2,errcts2,ext1,errext1,ext2,errext2,ml1,ml2,bkg1,bkg2,idsrc1,idsrc2\n') for var in variableSources: f.write(var+'\n') pass pass print('Wrote %s variable sources' %(len(variableSources))) fileroot = myslice.split('.out')[0] #header = pyfits.getheader('%s.img' %(fileroot),0) #filtered = header['XPROC0'].split('=')[1].split(':')[0] for cvar,var in enumerate(variableSources): outroot = '%s_var%s' %(var.split('.out')[0],cvar+1) boutroot = '%s_var%s_bkg' %(var.split('.out')[0],cvar+1) elements = var.split() ra,dec = (elements[1],elements[2]) selectEvents(filtered,ra,dec,20.0,outroot,20.0,30.0,boutroot) for ctrans,trans in enumerate(transients): outroot = '%s_trans%s' %(trans.split('.out')[0],ctrans+1) boutroot = '%s_trans%s_bkg' %(trans.split('.out')[0],ctrans+1) elements = trans.split() ra,dec = (elements[1],elements[2]) selectEvents(filtered,ra,dec,20.0,outroot,20.0,30.0,boutroot)