dir='/data/binary/summary/' ;-- lvl of confidence -> z* value ; lvl=80 & z=1.28 ; lvl=90 & z=1.645 lvl=95 & z=1.96 ; lvl=98 & z=2.33 ; lvl=99 & z=2.58 ;------------------------------------------------------------------------------------------ ;------------------------------------------------------------------------------------------ ;------------------------------------------------------------------------------------------ ;--- back readcol, dir+'taxo-pop.csv', class, bgNEA, bgMC, bgHun, bgMB, delimiter=',', /Silent, format='(A,F,F,F,F)' totBgNEA = total(bgNEA) totBgMC = total(bgMC ) totBgHun = total(bgHun) totBgMB = total(bgMB ) bgNEA /= totBgNEA bgMC /= totBgMC bgHun /= totBgHun bgMB /= totBgMB uBgNEA = z*sqrt( bgNEA*(1.-bgNEA) / totBgNEA ) uBgMC = z*sqrt( bgMC *(1.-bgMC) / totBgMC ) uBgHun = z*sqrt( bgHun*(1.-bgHun) / totBgHun ) uBgMB = z*sqrt( bgMB *(1.-bgMB) / totBgMB ) noNEA=where(bgNEA eq 0) if noNEA[0] ne -1 then uBgNEA[noNEA]=sqrt(totBgNEA)/totBgNEA noMC=where(bgMC eq 0) if noMC[0] ne -1 then uBgMC[noMC]=sqrt(totBgMC)/totBgMC noHun=where(bgHun eq 0) if noHun[0] ne -1 then uBgHun[noHun]=sqrt(totBgHun)/totBgHun noMB=where(bgMB eq 0) if noMB[0] ne -1 then uBgMB[noMB]=sqrt(totBgMB)/totBgMB print, '' print, ' BACKGROUND' print, 'Class NEA MC Hun MB' forprint, class, $ bgNEA*100, uBgNEA*100, $ bgMC*100 , uBgMC*100, $ bgHun*100, uBgHun*100, $ bgMB*100 , uBgMB*100, $ format='(A-4,4(4x,F4.1,2x,F4.1))' print, '#', totBgNEA, totBgMC, totBgHun, totBgMB, $ format='(A-4,4(7x,I4,3x))' print, '' ;------------------------------------------------------------------------------------------ ;------------------------------------------------------------------------------------------ ;------------------------------------------------------------------------------------------ ;--- Bin readcol, dir+'summary-bin.csv', class, binNEA, binMC, binHun, binMB, delimiter=',', /Silent, format='(A,F,F,F,F)' ; binNEA *= 10 ; binMB *= 10 totBinNEA = total(binNEA) totBinMC = total(binMC ) totBinHun = total(binHun) totBinMB = total(binMB ) ; print, 'Total binaries: ', totBinNEA, totBinMC, totBinHun, totBinMB binNEA /= totBinNEA binMC /= totBinMC binHun /= totBinHun binMB /= totBinMB uBinNEA = z*sqrt( binNEA*(1.-binNEA) / totBinNEA ) uBinMC = z*sqrt( binMC *(1.-binMC) / totBinMC ) uBinHun = z*sqrt( binHun*(1.-binHun) / totBinHun ) uBinMB = z*sqrt( binMB *(1.-binMB) / totBinMB ) noNEA=where(binNEA eq 0) if noNEA[0] ne -1 then uBinNEA[noNEA]=sqrt(totBinNEA)/totBinNEA noMC=where(binMC eq 0) if noMC[0] ne -1 then uBinMC[noMC]=sqrt(totBinMC)/totBinMC noHun=where(binHun eq 0) if noHun[0] ne -1 then uBinHun[noHun]=sqrt(totBinHun)/totBinHun noMB=where(binMB eq 0) if noMB[0] ne -1 then uBinMB[noMB]=sqrt(totBinMB)/totBinMB nbClass=n_elements(class) print, '' print, ' BINARY' print, 'Class NEA MC Hun MB' forprint, class, $ binNEA*100, uBinNEA*100, $ binMC*100 , uBinMC*100, $ binHun*100, uBinHun*100, $ binMB*100 , uBinMB*100, $ format='(A-4,4(4x,F4.1,2x,F4.1))' print, '#', totBinNEA, totBinMC, totBinHun, totBinMB, $ format='(A-4,4(7x,I4,3x))' print, '' relBinNEA = uBinNEA/BinNEA relBinMC = uBinMC /BinMC relBinHun = uBinHun/BinHun relBinMB = uBinMB /BinMB relBgNEA = uBgNEA/BgNEA relBgMC = uBgMC /BgMC relBgHun = uBgHun/BgHun relBgMB = uBgMB /BgMB print, '' print, ' RATIO' print, 'Class NEA MC Hun MB' forprint, class, $ binNEA/bgNEA, (binNEA/bgNEA)*sqrt( (relBinNEA)^2. + (relBgNEA)^2. ), $ binMC /bgMC , (binMC /bgMC )*sqrt( (relBinMC )^2. + (relBgMC )^2. ), $ binHun/bgHun, (binHun/bgHun)*sqrt( (relBinHun)^2. + (relBgHun)^2. ), $ binMB /bgMB , (binMB /bgMB )*sqrt( (relBinMB )^2. + (relBgMB )^2. ), $ format='(A-4,4(4x,F4.1,2x,F4.1))' print, '#', totBinNEA, totBinMC, totBinHun, totBinMB, $ format='(A-4,4(7x,I4,3x))' print, '' ;------------------------------------------------------------------------------------------ ;------------------------------------------------------------------------------------------ ;------------------------------------------------------------------------------------------ ;- Plot BINARY vs BACKGROUND yRang=[0,.6] cgPS_open, xSize=16, ySize=21, Filename=dir+'taxo-comparison.eps', $ /Metric, /Decomposed, /Landscape, /Encapsulated, Language_Level=2, /Quiet ticks=[' ', class, ' '] x=indgen(10)+1 xs=0.02 yRange=[0,65] xArr = indgen(nbClass) xShi=0.05 xBot = 0.15 yBot = 0.05 xLen = 0.84 yLen = 0.45 binNEA*=100 uBinNEA*=100 bgNEA*=100 uBgNEA*=100 binMB*=100 uBinMB*=100 bgMB*=100 uBgMB*=100 cgText, 0.05,0.5, /normal, orient=90, align=0.5, 'Fraction of the population (%)' botBinNEA = uBinNEA topBinNEA = uBinNEA botBgNEA = uBgNEA topBgNEA = uBgNEA low = where (binNEA-uBinNEA le 0) if low[0] ne -1 then botBinNEA[low] = BinNEA[low] hig = where (binNEA+uBinNEA ge yRange[1]) if hig[0] ne -1 then topBinNEA[hig] = yRange[1]-binNEA[hig] low = where (bgNEA-uBgNEA le 0) if low[0] ne -1 then botBgNEA[low] = BgNEA[low] hig = where (bgNEA+uBgNEA ge yRange[1]) if hig[0] ne -1 then topBgNEA[hig] = yRange[1]-bgNEA[hig] cgPlot, x+xs, binNEA, err_yLow=botBinNEA, err_yHigh=uBinNEA, psym='Open Circle', $ xTickName=replicate(' ',15), xtickint=1, xRange=[0,11], xminor=1, $ yRange=yRange, color='red', $ yminor=5, $ /Normal, /NoErase, $ position=[xBot, yBot+yLen, xBot+xLen, yBot+2*yLen] cgPlot, x-xs, bgNEA, err_yLow=botBgNEA,err_yHigh=topBgNEA, /OverPlot, psym='Open Square' xPop = 0.5 yPop = 58 sPop = 1.5 cgText, /Data, xPop, yPop, charSize=sPop, 'Near_Earth Asteroids' xSym = 1 xTxt = 1.5 sTop = 52 sBot = 44 sErr = 3 cgPlot, /OverPlot, xSym, sTop, err_yLow=sErr, err_yHigh=sErr, psym='Open Circle', color='Red' cgPlot, /OverPlot, xSym, sBot, err_yLow=sErr, err_yHigh=sErr, psym='Open Square', color='Black' cgText, /Data, xTxt, sTop-1.5, 'Binary' , charSize=sPop cgText, /Data, xTxt, sBot-1.5, 'Background', charSize=sPop cgAxis, xAxis=1, xRange=[0,11], xtickint=1, xtickname=ticks, xminor=1 botBinMB = uBinMB topBinMB = uBinMB botBgMB = uBgMB topBgMB = uBgMB low = where (binMB-uBinMB le 0) if low[0] ne -1 then botBinMB[low] = BinMB[low] hig = where (binMB+uBinMB ge yRange[1]) if hig[0] ne -1 then topBinMB[hig] = yRange[1]-binMB[hig] low = where (bgMB-uBgMB le 0) if low[0] ne -1 then botBgMB[low] = BgMB[low] hig = where (bgMB+uBgMB ge yRange[1]) if hig[0] ne -1 then topBgMB[hig] = yRange[1]-bgMB[hig] cgPlot, x+xs, binMB, err_yLow=botBinMB, err_yHigh=topBinMB, psym='Filled Circle', $ xTickName=ticks, xtickint=1, xRange=[0,11], xminor=1, $ yRange=yRange, color='red', $ yminor=5, $ /Normal, /NoErase, $ position=[xBot, yBot, xBot+xLen, yBot+yLen] cgPlot, x-xs, bgMB, err_yLow=botBgMB,err_yHigh=topBgMB, /OverPlot, psym='Filled Square' cgText, /Data, xPop, yPop, charSize=sPop, 'Inner Main Belt' cgPlot, /OverPlot, xSym, sTop, err_yLow=sErr, err_yHigh=sErr, psym='Filled Circle', color='Red' cgPlot, /OverPlot, xSym, sBot, err_yLow=sErr, err_yHigh=sErr, psym='Filled Square', color='Black' cgText, /Data, xTxt, sTop-1.5, 'Binary' , charSize=sPop cgText, /Data, xTxt, sBot-1.5, 'Background', charSize=sPop ; print, total(binNEA), total(binMB) cgPS_close, /png, Delete_PS=0 stop ;------------------------------------------------------------------------------------ ;------------------------------------------------------------------------------------ ;------------------------------------------------------------------------------------ ;------------------------------------------------------------------------------------ ;forprint, uBgNEA*100, uBinNEA*100 ;stop ; cgPlot, x+xs, binNEA, err_yLow=uBinNEA, err_yHigh=uBinNEA, psym='Filled Circle', $ xTickName=ticks, xtickint=1, xRange=[0,11], xminor=1, $ yRange=yRange, color='red' cgPlot, x-xs, bgNEA, err_yLow=uBgNEA,err_yHigh=uBgNEA, /OverPlot, psym='Filled Circle' read, trash, prompt='go to MC?' cgPlot, x+xs, binMC, err_yLow=uBinMC, err_yHigh=uBinMC, psym='Filled Circle', $ xTickName=ticks, xtickint=1, xRange=[0,11], xminor=1, $ yRange=yRange, color='red' cgPlot, x-xs, bgMC, err_yLow=uBgMC,err_yHigh=uBgMC, /OverPlot, psym='Filled Circle' read, trash, prompt='go to Hun?' cgPlot, x+xs, binHun, err_yLow=uBinHun, err_yHigh=uBinHun, psym='Filled Circle', $ xTickName=ticks, xtickint=1, xRange=[0,11], xminor=1, $ yRange=yRange, color='red' cgPlot, x-xs, bgHun, err_yLow=uBgHun,err_yHigh=uBgHun, /OverPlot, psym='Filled Circle' read, trash, prompt='go to MB?' cgPlot, x+xs, binMB, err_yLow=uBinMB, err_yHigh=uBinMB, psym='Filled Circle', $ xTickName=ticks, xtickint=1, xRange=[0,11], xminor=1, $ yRange=yRange, color='red' cgPlot, x-xs, bgMB, err_yLow=uBgMB,err_yHigh=uBgMB, /OverPlot, psym='Filled Circle' ; ; ;; read, trash, prompt='go to Summary?' ; ; ; yRange=[0,10] ; ; cgPlot, x-2*xs, binNEA/bgNEA, $ ; err_yLow= sqrt(uBinNEA^2 + uBgNEA^2), $ ; err_yHigh=sqrt(uBinNEA^2 + uBgNEA^2), psym='Filled Circle', $ ; xTickName=ticks, xtickint=1, xRange=[0,11], xminor=1, $ ; yRange=yRange, color='red' ; ; cgPlot, x-1*xs, binMC/bgMC, /OverPlot, $ ; err_yLow=sqrt(uBinMC^2 + uBgMC^2), $ ; err_yHigh=sqrt(uBinMC^2 + uBgMC^2), psym='Filled Square', color='Navy' ; ; cgPlot, x+1*xs, binHun/bgHun, /OverPlot, $ ; err_yLow=sqrt(uBinHun^2 + uBgHun^2), $ ; err_yHigh=sqrt(uBinHun^2 + uBgHun^2), psym='Filled Square', color='Green' ; ; cgPlot, x+2*xs, binMB/bgMB, /OverPlot, $ ; err_yLow=sqrt(uBinMB^2 + uBgMB^2), $ ; err_yHigh=sqrt(uBinMB^2 + uBgMB^2), psym='Filled Square', color='CornflowerBlue' ; ;stop ; ;------------------------------------------------------------------------------------------ ;------------------------------------------------------------------------------------------ ;------------------------------------------------------------------------------------------ ;- Ration binaries to bacrkground fracNEA = (1.*binNEA) / (bgNEA) fracMC = (1.*binMC ) / (bgMC ) fracHun = (1.*binHun) / (bgHun) fracMB = (1.*binMB ) / (bgMB ) uncNEA = fracNEA * sqrt( (uBgNEA/bgNEA)^2 + (uBinNEA/binNEA)^2 ) uncMC = fracMC * sqrt( (uBgMC/bgMC) ^2 + (uBinMC/binMC) ^2 ) uncHun = fracHun * sqrt( (uBgHun/bgHun)^2 + (uBinHun/binHun)^2 ) uncMB = fracMB * sqrt( (uBgMB/bgMB) ^2 + (uBinMB/binMB) ^2 ) ;- uncNEA = abs( uBinNEA/bgNEA) + abs(binNEA*uBgNEA/ ( bgNEA^2 ) ) ;- uncMC = abs( uBinMC /bgMC ) + abs(binMC *uBgMC / ( bgMC ^2 ) ) ;- uncHun = abs( uBinHun/bgHun) + abs(binHun*uBgHun/ ( bgHun^2 ) ) ;- uncMB = abs( uBinMB /bgMB ) + abs(binMB *uBgMB / ( bgMB ^2 ) ) binAst = binNEA+binMC+binMB bgAst = bgNEA+bgMC+bgMB uBinAst = uBinNEA+uBinMC+uBinMB uBgAst = uBgNEA +uBgMC +uBgMB fracAst = (1.*binAst) / (bgAst) uncAst = abs( uBinAST/bgAST) + abs(binAST*uBgAST/ ( bgAST^2 ) ) yRang=[0,6] cgPS_open, xSize=21, ySize=16, Filename=dir+'taxo-distro.eps', $ /Metric, /Decomposed, /Landscape, /Encapsulated, Language_Level=2, /Quiet cgPlot, indgen(nbClass), fltarr(nbClass)+1, color='light gray', linestyle=2, $ xTickName=class, xtickint=1, xminor=1, $ /Normal, $ xRange = [-.5,9.5], $ yRange = yRang, yMinor=4, $ position=[0.1, 0.1, 0.95, 0.95], $ yTitle='Ratio of binary to background' xArr = indgen(nbClass) xShi=0.05 col=['Black', 'Orange','Dark Gray','Cornflower blue'] sym=['FilledUpTriangle', 'Filled Circle', 'Filled Square', 'FilledDownTriangle'] shi=[-2,-1,0,1]*0.15 popName=['NEA','MC','Hungaria','IMB'] xKey = 0.5 yKey = 5.25 yShi = -0.5 yLen = 0.1 for kPop=0, 3 do begin case kPop of 0: begin x = xArr y = fracNEA u = uncNEA end 1: begin x = xArr y = fracMC u = uncMC end 2: begin x = xArr y = fracHun u = uncHun end 3: begin x = xArr y = fracMB u = uncMB end ; else: begin ; x = xArr ; y = fracMB ; u = uncMB ; end endcase ;--Check where values bad=where( ~finite(y) ) if bad[0] ne -1 then y[bad]=0 ;-Adjust uncertainties to graph low = u bad = where( (y-low) lt 0 ) if bad[0] ne -1 then low[bad]=y[bad] high = u bad = where( (y+high) gt yRang[1]) if bad[0] ne -1 then high[bad]=yRang[1]-y[bad] ;-Data points ok = where( y ge yRang[0] and y le yRang[1]) cgPlot, /OverPlot, x[ok]+shi[kPop], y[ok], psym=sym[kPop], $ err_yLow=low[ok], err_yHigh=high[ok], color=col[kPop] ;---Caption cgPlot, /OverPlot, xKey, yKey+kPop*yShi, psym=sym[kPop], $ err_yLow=yLen, err_yHigh=yLen, color=col[kPop] cgText, xKey+0.3, yKey+kPop*yShi-0.1, popName[kPop] endfor ; ;; nonZero = where( fracNEA ne 0 and bgNEA ne 0 ) ; nonZero = where( bgNEA ne 0 ) ; cgPlot, /OverPlot, xArr[nonZero]-2*xShi, fracNEA[nonZero], symSize=1.5, psym='Filled Circle', $ ; err_yLow=uncNEA[nonZero], err_yHigh=uncNEA[nonZero] ; ;; nonZero = where( fracMC ne 0 and bgMC ne 0) ; nonZero = where( bgMC ne 0) ; cgPlot, /OverPlot, xArr[nonZero]-1*xShi, fracMC[nonZero], symSize=1.5, psym='Filled Square', color='Orange', $ ; err_yLow=uncMC[nonZero], err_yHigh=uncMC[nonZero] ; ;; nonZero = where( fracHun ne 0 and bgHun ne 0 ) ; nonZero = where( bgHun ne 0 ) ; cgPlot, /OverPlot, xArr[nonZero]+1*xShi, fracHun[nonZero], symSize=1.5, psym='Filled Diamond', color='Gray', $ ; err_yLow=uncHun[nonZero], err_yHigh=uncHun[nonZero] ; ;; nonZero = where( fracMB ne 0 and bgMB ne 0 ) ; nonZero = where( bgMB ne 0 ) ; cgPlot, /OverPlot, xArr[nonZero]+2*xShi, fracMB[nonZero], symSize=1.5, psym='Filled Circle', color='Cornflower Blue', $ ; err_yLow=uncMB[nonZero], err_yHigh=uncMB[nonZero] ; ; nonZero = where( fracAst ne 0 ) ; cgPlot, /OverPlot, xArr[nonZero]+3*xShi, fracAst[nonZero], symSize=1.5, psym='Filled Circle', color='Red', $ ; err_yLow=uncAst[nonZero], err_yHigh=uncAst[nonZero] ; ; ; forprint, class, (1-fracNEA)/uncNEA, (1-fracMC)/uncMC, (1-fracHun)/uncHun, (1-fracMB)/uncMB ; forprint, class, fracMB, uncMB, uBgMB, uBinMB, bgMB ; print, 'class pop(%) unc bin(%) unc frac unc' ; forprint, class, bgMB, uBgMB, binMB, uBinMB, fracMB, uncMB; cgPS_close, /png, Delete_PS=0 print, 'NEA Q/S' QS = fracNEA[6]/fracNEA[7] uQS= sqrt( (uncNEA[6]/fracNEA[6])^2 + (uncNEA[7]/fracNEA[7])^2 ) print, QS, uQS, abs(1.-QS)/uQS print, '' print, 'NEA + MC + MB Q/S' QS = fracAST[6]/fracAST[7] uQS= sqrt( (uncAST[6]/fracAST[6])^2 + (uncAST[7]/fracAST[7])^2 ) print, QS, uQS, abs(1.-QS)/uQS print, '' print, 'NEA + MC + MB C/S' CS = fracAST[2]/fracAST[7] uCS= sqrt( (uncAST[2]/fracAST[2])^2 + (uncAST[7]/fracAST[7])^2 ) print, CS, uCS, abs(1.-CS)/uCS end