例によってこの文章を書いている今の時点では、寄る年波には勝てず記憶はきわめて怪しげなので、以下の(これまでも)記述には大嘘が多分に含まれるものと思われます。
DCCT装置からは、ハイインピーダンスのビーム電流信号が出てきています。1V=100mAではなかったかしら(あやしげ)。これをアドバンテスト社デジタルマルチメータR6551(DVM)で読みとり、GP-IBを通してEPICSに渡します。これを直接表示するとともに、シーケンサーでビーム寿命及び入射レートを計算します。
Device(R6551) %{ }% NameTable Functions = "DCV", "ACV", "2WR", "4WR", "DCA", "ACA"; NameTable DCVRange = "AUTO", "300mV", "3000mV", "30V", "300V", "1000V"; NameTable ACVRange = "AUTO", "300mV", "3000mV", "30V", "300V", "700V"; NameTable ResisRange = "AUTO", "300 Ohm"," 3000 Ohm", "30 kOhm", "300 kOhm", "3000 kOhm", "30 MOhm", "300 MOhm"; NameTable CurrentRange = "AUTO", "300mA", "3000mA"; NameTable SampleMode = "Free", "Hold"; NameTable SampleRate = "Fast", "Middle", "Slow"; NameTable ViewCols = "3+1/2", "4+1/2", "5+1/2"; NameTable Nulls = "OFF", "ON"; NameTable Scaling = "OFF", "ON"; NameTable Filter = "OFF", "ON"; NameTable AutoZero = "OFF", "ON", "ONCE"; EfastTable FunctionsSet = "F1","F2","F3","F4","F5","F6"; EfastTable DCVRangeSet = "F1,R0","F1,R3","F1,R4","F1,R5","F1,R6","F1,R7"; EfastTable ACVRangeSet = "F2,R0","F2,R3","F2,R4","F2,R5","F2,R6","F2,R7"; EfastTable TwResisRangeSet = "F3,R0","F3,R3","F3,R4","F3,R5","F3,R6","F3,R7","F3,R8","F3,R9"; EfastTable FwResisRangeSet = "F4,R0","F4,R3","F4,R4","F4,R5","F4,R6","F4,R7","F4,R8","F4,R9"; EfastTable DCCurrentRangeSet = "F5,R0", "F5,R6", "F5,R7"; EfastTable ACCurrentRangeSet = "F6,R0", "F6,R6", "F6,R7"; EfastTable SampleModeSet = "M0", "M1"; EfastTable SampleRateSet = "PR1","PR2","PR3"; EfastTable ViewColsSet = "RE3","RE4","RE5"; EfastTable NullsSet = "NL0","NL1"; EfastTable ScalingSet = "SC0","SC1"; EfastTable FilterSet = "FL0","FL1"; EfastTable AutoZeroSet = "AZ0","AZ1","AZ2"; ParamTable{ set_function {rec=mbbo, efast=FunctionsSet, name=Functions} set_DCV_range {rec=mbbo, efast=DCVRangeSet, name=DCVRange} set_ACV_range {rec=mbbo, efast=ACVRangeSet, name=ACVRange} set_TwR_range {rec=mbbo, efast=TwResisRangeSet, name=ResisRange} set_Fwr_range {rec=mbbo, efast=FwResisRangeSet, name=ResisRange} set_DCC_range {rec=mbbo, efast=DCCurrentRangeSet, name=CurrentRange} set_ACC_range {rec=mbbo, efast=ACCurrentRangeSet, name=CurrentRange} set_Samp_mode {rec=bo, efast=SampleModeSet, name=SampleMode} set_Samp_rate {rec=mbbo, efast=SampleRateSet, name=SampleRate} set_VCols {rec=mbbo, efast=ViewColsSet, name=ViewCols} set_Nulls {rec=bo, efast=NullsSet, name=Nulls} set_Scaling {rec=bo, efast=ScalingSet, name=Scaling} set_Filter {rec=bo, efast=FilterSet, name=Filter} set_Autozero {rec=mbbo, efast=AutoZeroSet, name=AutoZero} get_data {rec=ai, command="E\r\n",conv="%lf\r\n",leng=511} } %{ }%内容については説明するまでも無いと思います。なお、DCしか使いませんので、上の関数のうち使うのはget_data、set_DCV_range、set_SAMP_rage位です。
このほかの、例えばfit_startとかは、次のビーム寿命等の計算の時に使います。
I=I0exp(-t/τ)
この「τ」が(色々な効果を合わせた)ビーム寿命となる。実際にはARの真空系はあまりパワフルとは言えず(控えめな言い方)、ビーム電流の大きいところと小さな所ではかなり真空度が違い、ビーム寿命は上の式からかなり逸脱する(exponentailでなく直線的になる)が、ある区間を限ってみれば、exponentialの式でビーム寿命を定義しても問題はない。
実際には、1秒刻みでサンプルした今から過去に向かっての電流値のデータをいくつか記録しておいて、それをexponentialの式で最小2乗近似ををしてビーム寿命を出すことになります。多くの点数を取ればそれだけ1点毎のノイズ等の影響は小さくなりますが、あまり寿命が短いときに延々と長くフィット点数を取れば、元々の式からの逸脱も大きくなるし、また寿命の急変等にも対応出来なくなります。そこで、「今」のlifeの値をもとに寿命計算のフィット点数を判断する必要があります。また、この点数の切り替えにより(計算が変わるので)計算したlifeに若干の差がでることもあり、下手をすると切り替わりが非常に頻繁に起こってしまいます。そこで、切り替えにはシュミットトリガ的な動作が必要になります。シーケンサーの全ソースコードはここをクリックして下さい。
program ardcct_l2 %{extern double log();}% %{extern double fabs(); extern double exp();}% float FB; float LIFE; float RATE; %{ float ib[500]; float s_xy,s_x,s_y,s_x2,a1,a2,a3,a4,a5,a6; float samp_time = 1.0; int l_10, l_60, l_240; int llflag1,llflag2,llflag3,llflag4; }%logとかabsとかexpはシーケンサーとしてexternalファンクションになりますので、「c言語の」externalだとして宣言してやる必要があります。また、チャンネルと直接お話をしない部分は普通のCの書き方をすればよいので、%{と}%でくくって書きます。
assign FB to "fb_ar_dcct:dcct"; assign LIFE to "fb_ar_dcct:life"; assign RATE to "fb_ar_dcct:rate"; assign FITSTART to "fb_ar_dcct:fitstart"; assign P_MON01 to "fb_ar_dcct:dcct"; monitor P_MON01;チャンネルへのassignment及びmonitor宣言をします。今回のmonitorはdcct電流値の計算終了です(DVMへのサンプルと同じ)。これは実は値が変わらないと動かないので、本当はちょっとまずい。
ss state_main1 { state init { when() { printf("init state\n"); ZERO = 0.0; MAXFIT = 500; LIFE = 0.001; FITSTART = MAXFIT -10; FITN = MAXFIT -FITSTART; MIN=0.2; MIN2=0.003; P_MON01P=0; samp_time = 1.0; llflag1 = 0; llflag2 = 0; llflag3 = 0; llflag4 = 1; l_10 = 0; l_60 = 0; l_240 = 0; for (i=0; i<200; i++) { ib[i]=0.01; } }state sta1 }initステートは立ち上がったとき1回だけ動くsteteです。次はsta1ステートへ行き、戻っては来ません。
state sta1 { when(P_MON01!=P_MON01P) { s_x = 0.0; s_x2 = 0.0; s_y = 0.0; s_xy = 0.0; pvGet(FB); /* printf("FB= %f\n",FB); */ if (FB<MIN) ib[MAXFIT-1]=0.01; else ib[MAXFIT-1] = FB;最小2乗法で簡単に計算するために、ビーム電流のlogをとります。このため、電流値が0とか負の値は大変困ってしまうので、その時はあきらめて0.01mAに設定してしまいます。
if ((LIFE > 30.0) && (llflag3)) { if ((l_240 > 239) && (!llflag4)) { FITSTART=MAXFIT-240; pvPut(FITSTART); llflag4=1; l_240 = 240; } else { l_240++; } } else if (LIFE > 15.0) { if ((l_60 > 59) && (!llflag3) || (llflag4)) { FITSTART=MAXFIT-60; pvPut(FITSTART); llflag3 = 1; llflag4 = 0; l_60 = 60; l_240 = 60; } else { l_60++; } } else if (LIFE> 1.0) { llflag3 = 0; llflag4 = 0; l_60 = 0; l_240 = 0; FITSTART=MAXFIT-10; pvPut(FITSTART); } else if (LIFE> 0.0) { llflag3 = 0; llflag4 = 0; FITSTART = MAXFIT - 5; pvPut(FITSTART); l_60 = 0; l_240 = 0; } else { llflag3 = 0; llflag4 = 0; l_60 = 0; l_240 = 0; FITSTART = MAXFIT - 3; pvPut(FITSTART); }今のビーム寿命をもとにフィットする点数を決めます。基本的に、寿命が30分以上なら240点(4分)、15分以上なら60点(1分)、1分以上なら10点、それ以下でビームがふえていなければ5点、ビームがふえているなら(入射中のはず)3点としました。また、一気にモード変更が起きると拙いので、その寿命が何回か続くと次のモードへ移行するようにしました。ただし、lifeが急に長くなる理由はあまりないと思いますが、不幸にしてlifeが急に短くなる理由はいくらでもありますので、短くなったときは問答無用で短寿命モードへ移行するようにしました。今回のARの運転では、EPICSデータベース上では今がどんな運転モードか(入射中か、2.5GeV蓄積中か、加速中か、6.5GeV蓄積中か)分からないので、やはり入射近辺ではライフ計算モード移行がうまくいかず、しばらくがたがたが続くこともありました。この辺はなかなか奥が深そうで、下手をするとビーム寿命計算がライフワークになってしまうかも知れません(そいつはご勘弁を)。
for (i=FITSTART;i<MAXFIT; i++) { s_x += i*samp_time; s_x2 += (i*samp_time)*(i*samp_time); s_y += log(ib[i]); s_xy += (i*samp_time)*log(ib[i]); } FITN=MAXFIT-FITSTART; a1 = FITN*s_x2 - s_x*s_x; a2 = FITN*s_xy - s_x*s_y; a3 = a2/a1; RATE = a3*ib[FITSTART]; pvPut(RATE); a3 = -a3*60; /* printf(" a3 and ABS =%f %f\n",a3,fabs(a3)); */ if (fabs(a3)<MIN2) a4 = 0; else a4 =1.0/a3; LIFE = a4; pvPut(LIFE); P_MON01P=P_MON01; for (i=0;i<MAXFIT-1; i++) { ib[i] = ib[i+1]; } }state sta1 } }実際の最小2乗計算をしているところです。線形フィットをして、傾き分の1が寿命になります。これを分の単位に換算して終わりです。
実際には、運転中にシーケンサを止め、改造し、再ロードしたくなることが多くあると思われますが、不思議なことに今まで再ロードに成功した試しがありません。よくいって動作しないか、あるいはIOCがハングアップしてしまうといった状態です。まあ、あきらめてリブートする方が早いかも知れませんが、運転中にDCCTの情報が無くなるといった事態が起きて良い物か、どうでしょう。