memo391 Forward <= Today => Previous memo386

ふと見るとindigoのメインスクリプトのファイルサイズが65523バイトだった(2:40:44)

*WWW::Mechanizeであそぼう（端末からindigoに画像アップロード） => memo389
自分で作っておいて何だが久々のヒットの予感

*dTの端を検知する => memo389
*宇宙線ゲートを作る => memo389
*のために差分スペクトルを描こう => memo389

*indigoマルチユーザ化したい
*簡易encrypt-decrypterつけたい（これはすぐできそう？）
*kumacメインソースコード [cernlibを展開したもの]/src/packlib/kuip/code_kuip/kmacro.c

*端末からtwitterに投稿できるようにしたらもっとtwitter使うようになるかも。

*眼鏡が広角レンズのかわりになるらしい、と聞いたのでやってみたんだが、
微妙な差。

*'SELF'聴いてるときに椅子ずらしたらアーム削れる時の音がした(6:54:12)

<imagesender.pl>
リモートログイン時に一々手元にファイル転送してからアップロードするのが嫌なので、
nebula01のターミナルから直接indigoにファイルをアップロードする。
nebula01上ではコマンド is として使用できる。

*param
1st param -> image file path
2nd param -> title (空の場合はImagesenderが聞き返してくる)
3rd param -> comment (arb)
title, commentは日本語もだいたいいける。

*option
-n -> no comment
   コマンドラインでコメントを入れなくても聞き返されない
-v -> verbose mode
   ことあるごとに喋るモード

*usage
アップするファイル名だけ指定すればとりあえず動く。
アップロード成功の場合は、作成されたクリップのタイトルが表示される。
失敗したときは失敗メッセージの後にHTTPステータスコードが表示される。

アップロードテスト


端末入力:
[nebula@nebula01 kawada]$ ./imagesender.pl imgs/horns.ps -v
imgs/horns.ps
Enter title: uptest2
Any comment: uppppppppppppppp
Sending file 'imgs/horns.ps'...
Caption is 'uptest2'
Full path: '/home/nebula/exp/test-20100924/users/kawada/imgs/horns.ps'
Uploaded Successfully!!
Indigo: clip389 uptest2
[nebula@nebula01 kawada]$

<centroid.kumac>
-------------------------------------------------------

 macro centroid thresh=100 bin=400 min=1 max=400 id inf=-40. sup=40.
  v/del par
  v/cre hist([bin]) R [bin]*0
  v/cre par(2) R
  hi/get_vect/content [id] hist
  j=0
  
  zone 1 2
  
  accum = 0
  chs   = 0
  do i=[min],[max]
    v = hist([i])
    if [v]>[thresh] then
      accum = [accum] + hist([i])
      chs   = [chs] + 1
    endif
  enddo
  
  avrg = [accum]/[chs]
  avrg = [avrg]/2
  mess half active-average: [avrg]

  i = [min]
loop:
 v = hist([i])
 if [avrg]>[v] then
   i = [i] + 1
   goto loop
 endif
  mess exceed:[v] at [i]
 imin = [i]
 
loop2:
 v = hist([i])
 if [v]>[avrg] then
   i = [i] + 1
   goto loop2
 endif
 mess down:[v] at [i]
 imax = [i]
 
 v = [imin]*([sup]-[inf])/[bin]+[inf]
 vec/inp par(1) [v]
 mess [v]
 v = [imax]*([sup]-[inf])/[bin]+[inf]
 vec/inp par(2) [v]
 mess [v]
 
 
 
 v/del hist
return

macro many thresh=1000 bin=400 min=1  max=400 from num=1 inf=-40. sup=40.
  vec/cre rising([num]) R
  vec/cre falling([num]) R
  vec/cre ids([num]) R
  do i=1,[num]
    id=[from]+[i]-1
    exec kumac/centroid#centroid [thresh] [bin] [min] [max] [id] [inf] [sup]
    vec/inp ids([i]) [id]
    vec/inp rising([i]) par(1)
    vec/inp falling([i]) par(2)
    wait
  enddo
  vec/write ids,rising,falling ! '(f4, 1x, f8.2, 1x, f8.2)'
*  vec/write storepeak ./ycentroid.txt '(f8.2)'
return

<ゲートを作る>

宇宙線ピーク部分だけをオートマチックに抜きたいのだが。。。


この図でいうところの300-900あたりの山を抜きたい。
条件としては、
「下から見ていって、スペクトルがフラットな部分から急速に上がるところをマーク
->下り部分でそこでの高さと同じになるところをマーク」
２つのマーキングの間でゲートをかける。

とりあえず差分スペクトルを書かせる。

N101A*の差分スペクトル(上下でスケールが違う)

まじめな差分だとデジタルすぎてダメなので、
スムージングパラメータSを使って
f'(x) = {Σ{i=1,S}f(x+i)-f(x-i)}/S とする。

差分グラフで「ゼロ付近まで落ちてから急上昇するエリア」を選んでそのマイナスピークを抜けば良い。
スレッショルドによる崖を拾わないために、だいたい100ch以下は扱わないことにする。

********************************************************************
別のアプローチ
1) A*をAcalする(pedestal=0, 30MeV=peak)
2) 全neutでAcalのピークが同じだと思うと全neutで同じゲートが使える
3) コモンゲートでdTを抜く
4) posCal!

1)のためにA*のピークを抜く。
300-700で差分スペクトルが上からゼロクロスする点という感じでよかろう。
dq.kumacでピークを抜く。
それを使ってAcalするのはanapawの領分。

差分スペクトル描画マクロ
dq.kumac
-----------------------------------------------------------
 macro dq bin=4100 min=1 max=4100 id smooth=50 from=300 to=700
   v/del par
   v/cre xx([bin]) R
   sigma xx=array([bin],1#[bin])
   v/cre hist([bin]) R [bin]*0
   v/cre dhist([bin]) R [bin]*0
   v/cre par(2) R
   hi/get_vect/content [id] hist
   j=0
   zone 1 2
   hi/pl [id]
   max = [max]-[smooth]
   min = [min]+[smooth]
   
   v=0
   do j=1,[smooth]
     i2 = [min]+[j]
     i3 = [min]-[j]       
     v = [v] + hist([i2]) - hist([i3])
   enddo

   do i=[min],[max]
     v2 = [v] / [smooth]
     v/inp dhist([i]) [v2]
     i2 = [i]+1
     i3 = [i]-[smooth]
     v = [v] - hist([i2]) + hist([i3])
     i2 = [i]+[smooth]+1
     i3 = [i]
     v = [v] + hist([i2]) - hist([i3])
   enddo

   
   graph [bin] xx dhist
   
   do i=[from],[to]
     a = dhist([i])
     if [a]<0 then
       peak = [i]
       breakl
     endif
   enddo
   
   mess [peak]
 return

<Acal>

とりあえず


neut108を代表として dE*ゲートを作る。
---------------------------
ANAPAW> ht 108
 Current Histogram ID :  108
ANAPAW> xval; xval;
  
 ID:   108 | X Value:    15.98 | Contents:       2262
  
  
 ID:   108 | X Value:    50.97 | Contents:       2187
---------------------------
ゲートは[16:51]とする。

紆余曲折あってYキャリブレーション。

exper.kumac
減衰長抽出器（暫定版）
これでN101の減衰長を出したら3800とかだった。マジすか。
-------------------------------------------------------------------------
 macro exper bin1 min1 max1 bin2 min2 max2 id
   vec/cre hist2([bin1],[bin2]) R
   v/cre par(3) R
   his/get_vect/content [id] hist2
   vec/cre hist1([bin1]) R
   v/cre xx([bin1]) R

   v = hist2(60,60)
   mess [v]

   do i=1,[bin1]
     v = 0
     d = 0
     do j=1,[bin2]
       v = [v] + hist2([i],[j])*(([max2]-[min2])*[j]/[bin2] + [min2])
       v2 = hist2([i],[j])
       if hist2([i],[j])>0 then
* mess [v2]
         d = [d] + hist2([i],[j])
       endif
     enddo
     if [d]>0 then
       v = [v] / [d]
     else
       v = 0
     endif
     v/inp hist1([i]) [v]
     v2 = [v]
     v = [i]*([max1]-[min1])/[bin1]+[min1]
* mess [v] [v2]
     v/inp xx([i]) [v]
   enddo

* assuming neut is [-900,+900].
   downside = [bin1]*(-900-[min1])/([max1]-[min1])
   upside   = [bin1]*(+900-[min1])/([max1]-[min1])
   v = hist1([downside])
   v2= hist1(  [upside])

   bin3 = [upside]-[downside]+1

   v/cre hist3([bin3]) R
   v/cre xxx([bin3]) R
   do i=[downside],[upside]
     k = [i]-[downside]+1
     v = hist1([i]);
     v/inp hist3([k]) [v]
     v = [i]*([max1]-[min1])/[bin1]+[min1]
     v/inp xxx([k]) [v]
   enddo

   graph [bin3] xxx hist3
   vec/fit xxx hist3 ! e+p0 ! 3 par
   lambda = par(2)
   lambda = 1 / [lambda]
   mess attenuation length: [lambda]
 return