1

ROOT 解析入門 II∼ アホの子でもわかる物理解析基礎 ∼H. Ito2017 年 7 月 9 日

目次
1

はじめに

2

2

相対論入門

2

2.1

TGraph で関数をプロットしてみよう . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

22.2

粒子の運動量と速度の関係. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

2.3

飛行時間測定法 (Time of flight). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

52.4

粒子識別 (PID) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

粒子衝突と崩壊 ∼ 静止 K + 崩壊を中心に ∼

12

3.1

２体崩壊 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

3.2

3 体崩壊 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

3

3.2.1

K → 3π . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

3.2.2

K → πlν . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

3.3
4

静止 K 崩壊荷電粒子のスペクトラム . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

まとめ

36

2 相対論入門

1 はじめに
内容は素粒子実験物理学に絡めて、例題で ROOT 解析の手助けになればいいかなと思ってます。解析のコンセプ
トはスクリプトマクロ形式でコマンド root -l xxx.cxx だけで実効できるようにしています。重要なのはマクロの名前です。コードの中で void sample1() を定義していれば、sample1.cxx とファイル名にしないと ROOT から怒られますので注意してください。環境は ROOT ver. 5.34 を使用しています。2016 年ごろに ver.6 がリリースされましたが、未だにこちらを使っています。これ以降の話は C/C++ と Shell Script の内容が混じっている可能性がありま
すが、ちょいちょい説明し忘れることがあるので注意してください。インストールは別紙「アホの子でもできる n 次元解析」を参考にしてください。自分がアホでないと思っている方は是非閲覧しないでください。著者のアホ度がバレてしまうので。

2 相対論入門
粒子の運動学的な計算は学部生にとっては必修だと思います。
段階を経て理解できるように努めています。

2.1 TGraph で関数をプロットしてみよう
ここでは単純に数学の勉強です。サンプルを下に示します。

sample1.cxx

tiltley-axisvoid sample1(){
int n=100;

2
1.8
1.6

double x[100],y[100];

1.4

for(int i=0;i<100;i++){1.2
1x[i]=i*0.01;

0.8

y[i]=2*x[i];

0.6
0.4

}

0.2

graph=new TGraph(n,x,y);00graph->SetTitle("tiltle");

0.2

0.4

図1

0.6

0.8

1
x-axisy = 2x のグラフ

graph->GetXaxis()->SetTitle("x-axis");graph->GetYaxis()->SetTitle("y-axis");c1=new TCanvas("c1");graph->Draw("al");c1->Print("output.pdf");
return;}

これは y = 2x の曲線をプロットするプログラムです。100 回 for 文でループさせて xy の配列に格納させます。その後、TGraph を定義して、タイトル、x 軸 y 軸のラベルを導入します。TCanvas でキャンバスを定義してそこにgraph を描写する流れです。

2

2.2

2 相対論入門

粒子の運動量と速度の関係

2.2 粒子の運動量と速度の関係
静止系 K から速度比 β = v/c で動く系 K ′ へローレンツ変換した場合、粒子エネルギーと運動量の関係は

(

E′
p′

)

(
=

γ
γβ

γβ
γ

)(

)

E
p

(1)

で表せます。粒子静止系 K: p = (E, −p) = (m, 0) とすると、話がラクになりますね。ここから E = γm, p = γβm√
1 + (m/p)2 が導き出せます。ここで大事なのが粒子の質量です。

なので β = 1/2

• 電子質量 0.511 MeV/c ,2• µ 粒子質量 105.6 MeV/c ,2• π + 中間子質量 139.6 MeV/c ,2• K + 中間子の質量 493.7 MeV/c ,• 陽子質量 938.3 MeV/c

2を前提にしてサンプルソースを改良して運動量と速度比の関係をグラフにします。

sample2.cxx

void sample2(){
int n=100;
double x[100],y[100];
double m[5];m[0]=0.511;

ratio: β = v/c

Relation of momentum and velocity ratio
1 e
0.9

µ
π+

0.8

K

+

0.7
p

0.6
0.5

m[1]=105.6;

0.4m[2]=139.6;

0.3m[3]=493.7;

0.2m[4]=938.3;

0.1
0
0

200

400

600

800

1000

1200

1400

1600 1800 2000
momentum (MeV/c)graph[5]=new TGraph();for(int j=0;j<5;j++){for(int i=0;i<100;i++){x[i]=i*20+10;//p (MeV/c)y[i]=1./sqrt(1.+(m[j]/x[i])**2);//beta}
graph[j]=new TGraph(n,x,y);graph[j]->SetLineColor(j+2);
graph[j]->SetLineWidth(3);
}




3

2.2

2 相対論入門

粒子の運動量と速度の関係

sample2.cxx の続き

c1=new TCanvas("c1");
c1->SetGridx();c1->SetGridy();graph[0]->SetTitle("Relation of momentum and velocity ratio");graph[0]->GetXaxis()->SetTitle("momentum (MeV/c)");graph[0]->GetYaxis()->SetTitle("ratio: #beta = v/c");graph[0]->GetHistogram()->SetMaximum(1);graph[0]->GetHistogram()->SetMinimum(0);graph[0]->GetXaxis()->SetLimits(0,2000);graph[0]->Draw("al");for(int j=1;j<5;j++){graph[j]->Draw("l");}text=new TLatex(50,0.96,"e");text->SetTextColor(2);text->SetTextSize(0.03);
text->Draw("same");text=new TLatex(100,0.92,"#mu");text->SetTextColor(3);text->SetTextSize(0.03);
text->Draw("same");text=new TLatex(300,0.81,"#pi^{+}");
text->SetTextColor(4);text->SetTextSize(0.03);
text->Draw("same");text=new TLatex(500,0.75,"K^{+}");
text->SetTextColor(5);text->SetTextSize(0.03);
text->Draw("same");text=new TLatex(600,0.6,"p");text->SetTextColor(6);text->SetTextSize(0.03);
text->Draw("same");c1->Print("output.pdf");
return;
}




4

2.3

飛行時間測定法 (Time of flight)

2 相対論入門

ここで注目してもらいたいのは graph[5] が先に定義されていることです。本当は TGraph*graph[5] と定義するのでしょうが、ROOT さんはそこらへんが結構いい加減みたいです。すると graph[i] は i でループさせることができる
ので、粒子ごとにプロットができるようになります。そして、曲線の色は SetLineColor(int fColor) で引数 fColor の番号で色を指定します。TCanvas にグリッド線を入れてグラフが見やすくするのは趣味ですね。Log スケールにしたい場合は SetLogy() など挿入するといいでしょう。

TGraph の Draw() 関数の引数には”al”とありますが、”a”は軸を用意すると言う意味で新規でキャンバスに描写します。”l”はラインで折れ線グラフモード、他には”p”は点グラフになるので覚えておくと良いでしょう。
プロットするグラフの表示範囲を選択したい時に便利なコマンドを以下に示します。

graph->GetXaxis()->SetLimits(double min,double

max);graph->GetHistogram()->SetMaximxum(double max);
graph->GetHistogram()->SetMinimum(double min);
これはコマンドの名前通りで、X 座標のとりうる範囲を選択します。Y 軸は GetYaxis でとってくないで、GetHis-

togram でとってくるところが要注意です。
キャンバス内にテキストを表示させる関数は TText() もしくは TLatex() がある。著者のオススメは TLatex で

Tex の質感でギリシャ文字や特殊文字をキャンバスに代入できます。書き方の例は
text=new TLatex(100, 300, "test e #mu #nu #pi K p");text->SetTextSize(0.03);
text->SetTextColor(2);text->Darw("same");
です。TLatex の引数は x 座標、y 座標、そして文字の順番です。挿入する位置は文字 BOX をイメージすると左上の点に位置します。また、ポインタである text は本当なら、TLatex*text; とするところですし、何度も重複するのは良くないのですが、何度も言いますが、ROOT さんは適度にいい加減なので、何度同じポインタを使っても重ね書
きしてくれます。

2.3 飛行時間測定法 (Time of flight)
より実験的なグラフを書いてみましょう。粒子の速度が運動量に依存しているこということはわかりました。例えば、距離 L = 10m だけ離した 2 つの粒子検出器を用意し、片方をスタート、他方をストップで同じ運動量 p を持つ粒子を通過させます。運動量が同じで質量が違うので光速との時間差を観測することで粒子を特定することが可能です。この測定法は飛行時間測定法 (TOF: Time of flight) と呼ばれ粒子識別の一つとして現在の素粒子・原子核物理
学実験で採用されている技術です。

sample3.cxx

void sample3(){
int n=100;
double x[100],y[100];
double m[5]={0.511,105.6,139.6,493.7,938.3};
double c=2.99e8;
double L=30;//difference distance




5

2.3

飛行時間測定法 (Time of flight)

2 相対論入門sample3.cxx の続きRelation of momentum and time differencefor(int j=0;j<5;j++){for(int i=0;i<100;i++){x[i]=i*20+10;//p (MeV/c)time difference (ns)graph[5]=new TGraph();y[i]=(1.-

40
35
30
25

p

20

1./sqrt(1.+(m[j]/x[i])**2))

K

15/c*L*1e9;//dT(ns)}10graph[j]=new TGraph(n,x,y);5graph[j]->SetLineColor(j+2);

0
0

π

+

+

µ
e
200

400

600

800

1000

}

1200

1400

1600 1800 2000
momentum (MeV/c)

c1=new TCanvas("c1");
c1->SetGridx();c1->SetGridy();graph[0]->SetTitle("Relation of momentum and time difference");graph[0]->GetXaxis()->SetTitle("momentum (MeV/c)");graph[0]->GetYaxis()->SetTitle("time difference (ns)");graph[0]->GetHistogram()->SetMaximum(40);graph[0]->GetHistogram()->SetMinimum(-1);
graph[0]->GetXaxis()->SetLimits(0,2000);graph[0]->Draw("al");for(int j=1;j<5;j++){graph[j]->Draw("l");}text=new TLatex(50,0.5,"e");text->SetTextColor(2);text->SetTextSize(0.04);
text->Draw("same");text=new TLatex(200,5,"#mu");text->SetTextColor(3);text->SetTextSize(0.04);
text->Draw("same");text=new TLatex(400,10,"#pi^{+}");
text->SetTextColor(4);text->SetTextSize(0.04);
text->Draw("same");




6

2.4

粒子識別 (PID)

2 相対論入門sample3.cxx の続きtext=new TLatex(600,15,"K^{+}");
text->SetTextColor(5);text->SetTextSize(0.04);
text->Draw("same");text=new TLatex(1000,20,"p");text->SetTextColor(6);text->SetTextSize(0.04);
text->Draw("same");c1->Print("output.pdf");
return;}

出力結果の図を見ると、電子は数百 MeV でほぼ光速と等しく、粒子の質量が重たくなるにつれて時間差が大きくなります (当たり前でしょうが)。高エネルギー素粒子物理の領域では数 GeV/c の運動量領域で実験が行われることが多く、となると、時間差は電子と π 中間子で 1 ns も時間差が生まれないので、TOF で粒子識別することがほぼ無理。運動量が大きくなればなるほど TOF による粒子識別が困難になることがわかると思います。

2.4 粒子識別 (PID)
質量から粒子を同定することなどを粒子識別 (PID: particle identification) と呼びます。TOF について理解した
ところで、例題：1 GeV/c における TOF での π/K 識別はどのようにして行うのか？ まずデータがないのはいただけませんので、仮想的に生成しましょう。
運動量分解能が標準偏差 5 MeV/c で、時間分解能が標準偏差で 200 ps であると仮定します。100-1000 MeV/c の領域で π と K 中間子を発生させます。

sample4.cxx

PID

double mpi=139.6;
double mK=493.7;

dT (ns)void sample4(){14
12

double p;
double dt;

10

double c=3e8;//(m/s)

8double L=30;

64

his1=new TH2F("his1","his1",2

1e2,1e3,5e3,1e2,-1,100);
his2=new TH2F("his2","his2",0

1e2,1e3,5e3,1e2,-1,100);

1000



7

1500

2000

2500

3000

3500

4000

4500
5000
p (MeV/c)2.4

粒子識別 (PID)

2 相対論入門sample4.cxx の続きofstream ofs("data.dat");c1=new TCanvas("c1","",500,400);

for(int i=0;i<100000;i++){c1->SetGridx();

p = gRandom->Uniform(1000,5000);

c1->SetGridy();dt=(1.-1./sqrt(1.+(mpi/p)**2))/
c*L*1e9;his1->SetTitle("PID");

p=gRandom->Gaus(p,5);

his1->SetXTitle("p (MeV/c)");dt=gRandom->Gaus(dt,0.2);

his1->SetYTitle("dT (ns)");ofs<<p<<" "<<dt<<endl;

his1->GetXaxis()->SetRangeUser(1000,5000);

his1->Fill(p,dt);

his1->GetYaxis()->SetRangeUser(-1,15);
his1->Draw("box");dt=(1.-1./sqrt(1.+(mK/p)**2))/
c*L*1e9;his2->SetLineColor(2);p=gRandom->Gaus(p,5);his2->Draw("boxsame");dt=gRandom->Gaus(dt,0.2);
ofs<<p<<" "<<dt<<endl;

ofs.close();his2->Fill(p,dt);

c1->Print("output.pdf");}
return;gStyle->SetOptStat(0);

}





乱数は gRandom で生成することができます。

Uniform(start, stop) で start から stop までの領域を一様乱数で値を出力します。
例えば、

gRandom->Uniform(0,1);
では 0〜1 までの少数を一様に生成するのです。
一方、

gRandom->Gaus(mean, sigma)
は平均値 mean から標準偏差 sigma の正規分布の重みで乱数を生成します。

2 次元ヒストグラムの Draw オプションはいくつかあり、よく使用するのは box と colz, そして lego でしょう。重ね書きするときは same を加えると良いです。何もオプションをつけないと、ただの点としてプロットされます。そして紺色ヒストグラムが π で、赤ヒストグラムが K を示してます。

ofstream はファイル出力の C++ のクラスです。ofs() の引数に出力先を入力し、ofs<<a<<" "<<b<<endl;
で１行に a と b を出力することができます。ここで endl は改行を意味します。
データは π と K を交互に出力され、配列は運動量 (MeV/c) と時間差が格納されています。
ここで生成したデータを使って質量分布を作っていきましょう。

8

2.4

粒子識別 (PID)

2 相対論入門

system 関数は Linux 標準コマンドを実行してくれるので役立ちます。
wc コマンドは ファイルの行数を返すコマンドです。それにパイプをつないで、cat を命令します。cat は出力をそのまま返すのですが、> でつなぐとその後のファイルに出力出来ます。つまり、.lg というファイルに data.dat の行数
つまり、イベント数が格納されることになります。こうすることで、for 文でループを回す最大値を知ることができます。

sample5 に解析コードを示します。

sample5.cxx

void sample5(){
double m,m2;

2400
2200
2000
1800
1600
1400
1200
1000
800
600
400
200
0

double p;
double dt;
double c=3e8;//(m/s)

π+

−400

600

0

200

π+

1000

double L=30;

−200

400

K

600

800

char name[100];

−400

−200

0

200

π+

700
600

600

800

200

−400

ifstream lg(".lg");
400
300

0

200

π+

500

int loopmax=0;

−200

400

K

600

800

× 10

3

× 10

3

+

PID π /K: 1.57σ
PID K/ π : 1.38σ

1500 < p/(MeV/c) < 2000

200

lg>>loopmax;

3

+

100
0

× 10

PID π /K: 2.26σ
PID K/ π : 2.00σ1000 < p/(MeV/c) < 1500

300

system("wc -l data.dat|cat > .lg");

400

K

500
400

3

+

200
0

× 10

PID π /K: 3.60σ
PID K/ π : 3.22σ

500 < p/(MeV/c) < 1000

400

int region;+

PID π /K: 6.84σ
PID K/ π : 5.60σ

0 < p/(MeV/c) < 500

1200

800

K

100
0

−400

−200

0

200

400

600

800

lg.close();ifstream ifs("data.dat");
TH1S*mass[10];for(int j=0;j<8;j++){sprintf(name,"mass%d",j);mass[j]=new TH1S(name,name,2e2,-5e5,1e6);}
for(int i=0;i<loopmax;i++){ifs>>p>>dt;
beta=1.-dt*1e-9*c/L;
m2=p*p/beta/beta*(1+beta)*(1-beta);//MeV^2
region=(p-1e3)/5e2;
if(region>7) region =7;
mass[region]->Fill(m2);}
ifs.close();




9

2.4

粒子識別 (PID)

2 相対論入門sample5.cxx の続き

c1=new TCanvas("c1","",500,400);->SetLabelSize(0.06);

c1->Divide(1,4);gStyle->SetOptStat(0);sprintf(name,for(int j=0;j<4;j++){"%d < p/(MeV/c) < %d",c1->cd(j+1)->SetGridx();

j*500,(j+1)*500);

c1->cd(j+1)->SetGridy();text=new TLatex(-450e3,mass[j]->GetXaxis()mass[j]->GetMaximum()*0.5,name);->SetRangeUser(-5e5,8e5);

text->SetTextSize(0.12);

mass[j]->Draw();

text->SetTextColor(2);func=new TF1("func","gaus(0)+gaus(3)");text->Draw("same");func->SetLineStyle(2);func->SetParameters(

text=new TLatex(

2e3,1.5e4,5e4,2e3,2.4e5,5e4);-200e3,mass[j]->Fit("func","Q","",-1e5,4e5);mass[j]->GetMaximum()*0.8,"#pi^{+}");

func1=new TF1("func1","gaus",-5e5,1e6);

text->SetTextSize(0.2);

func1->SetLineColor(3);text->SetTextColor(3);func1->SetLineStyle(2);text->Draw("same");func1->SetParameters(
func->GetParameter(0),text=new TLatex(

func->GetParameter(1),400e3,

func->GetParameter(2)mass[j]->GetMaximum()*0.8,);"K^{+}");

func1->Draw("same");text->SetTextSize(0.2);

func2=new TF1("func2",

text->SetTextColor(4);"gaus",-5e5,1e6);

text->Draw("same");func2->SetLineColor(4);func2->SetLineStyle(2);double sigma=

func2->SetParameters((func2->GetParameter(1) -func->GetParameter(3),func1->GetParameter(1))/

func->GetParameter(4),func1->GetParameter(2)/2;

func->GetParameter(5)sprintf(name,);"PID #pi/K: %0.2f#sigma",sigma);func2->Draw("same");text=new TLatex(400e3,
mass[j]->GetMaximum()*0.5,name);mass[j]->SetTitle("");

text->SetTextSize(0.15);

mass[j]->GetXaxis()

text->SetTextColor(2);->SetLabelSize(0.1);

text->Draw("same");mass[j]->GetYaxis()



10

2.4

粒子識別 (PID)

2 相対論入門sample5.cxx の続きsigma =

func2->GetParameter(1)-(func2->GetParameter(1) -func1->GetParameter(1))/2,

func1->GetParameter(1))/

0,

func2->GetParameter(2)/2;(func2->GetParameter(1)-sprintf(name,func1->GetParameter(1))/2,"PID K/#pi: %0.2f#sigma",sigma);mass[j]->GetMaximum());

text=new TLatex(400e3,

line->SetLineColor(2);mass[j]->GetMaximum()*0.3,name);line->Draw("same");text->SetTextSize(0.15);}

text->SetTextColor(2);c1->Print("output.pdf");

text->Draw("same");return;}line=new TLine((

事前に mass という配列ポインタに TH1S のヒストグラムのクラスを付与しておきます。

TF1S*mass[10];その上で、for 文で new TH1S() でヒストグラムを定義します。
そして、ifstream を使ってデータを読み出します。m2 は質量の２乗を意味しています。mass は運動量ごとに Fill することで識別能力が運動量ごとにどう変わっていくか見ることができます。
フィッティングは２つの正規分布でフィットを行いますが、その前に２つのガウスを持つ関数を定義します。

func=new TF1(char name, char function);で定義します。一つ目の引数は関数名、２つ目は関数を char 型で書きます。ROOT の中には gaus という関数が事前に登録されていて、正規分布を示します。括弧の中身は自由パラメータが 0 から始まり、3 つ使用するので、gaus(0) = [0] ∗ exp(−((x − [1]/2/pi/[2])2 ))を示します。そして function の初期パラメータを設定します。最後にフィットをしますが、mass-¿Fit(char func-

tionname, char option1, char option2, double start, double stop); 最初の引数はフィットする関数名、option1 にQ を入れると、標準出力で情報を出さないようになります。基本的に option1, 2 はダブルクォーテーションでくくって、中身なしにするといいと思います。start と stop はフィット領域の最初と最後を示します。
粒子質量が正規分布で測定できたとして、ちょうど π と K のピークの半分を π/K 識別の境目にします。π の山のピークと境目が標準偏差の何倍かということで粒子識別能力を評価します。0-500 MeV/c の運動量領域では 6.84σ

500-1000 MeV/c の運動量領域では 3.60σ
1000-1500 MeV/c の運動量領域では 2.26σ
1500-2000 MeV/c の運動量領域では 1.57σ
この値が大きいと粒子識別能力が高いことを意味します。統計学でよく言われる標準偏差では 1 σが 68%,2 σが95%,3 σが 99.7% をカバーすることが知られています。これを指標に 3σ 以上あれば十分な識別能力があると判断す
るなら、TOF による粒子識別の限界は 1000 MeV/c あたりにあることがわかります。

11

3 粒子衝突と崩壊 ∼ 静止 K + 崩壊を中心に ∼3 粒子衝突と崩壊 ∼ 静止 K + 崩壊を中心に ∼おなじみの K + です。
これは us̄ で結合している中間子 (meson) ですよね。
質量は mK = 493.677 MeV.
寿命は τ = 12.38 ns.
主な崩壊チャンネルは表 1 に示します。
表1

Main K + decay channel

Decay channel

Branching ratio

K + → e+ νe (Ke2 )1.58 × 10−5

K + → µ+ νµ (Kµ2 )63.55%

K+→ e π νe (Ke3 )

5.07%

K

+

+ 0

→ µ π νµ (Kµ3 )

3.35%

K

+

→ π π (Kπ2 )20.66%

K

+

+ 0 0

→π π π

1.76%

K

+

+ + −

5.59%

+ 0

+ 0

→π π π3.1 ２体崩壊
ターゲット内に K + を止めて２対崩壊した場合、２粒子は Back to Back で放射するため、理解が簡単でしょう。
22
2Ke2 と Kµ2 はヘリシティ抑制によって前者が抑制されている典型的な例です。分岐幅は Γ = gl G
8π fK mK ml (1 −(ml /mK )2 )2 で与えられるため、同じセミレプトニック崩壊過程である Ke2 と Kµ2 の崩壊分岐比は質量の関係から計算できます。

RK =Γ(Ke2 )me 2=
Γ(Kµ2 )
mµ 2(

mK 2 − me 2
mK 2 − mµ 2)2

= 2.57 × 10−5
実際には制動放射による補正が入り、標準模型では RK SM = (2.477 ± 0.001) × 10−5 と計算されます。実験ではRK exp = (2.488 ± 0.009) × 10−5 が報告されています。
Kπ2 は単色の π 0 が生成されるため、π 0 → 2γ 事象から 2 つのγ線のエネルギーとそれぞれの運動量を観測するとピーク構造を得られます。そのため電磁カロリメータの較正測定によく使用されます。

2 体崩壊はエネルギーと運動量保存則、そして娘粒子２つの運動量絶対値が等しいことから、Ke2 , Kµ2 , Kπ2 のそれぞれ放射される荷電粒子のエネルギーと運動量が計算できます。

mK 2 = (El + Eν )2√√mK = p2 + m1 2 + p2 + m2 2
ここで、娘粒子の１つが ν なら質量 m2 = 0 と近似して運動量を決定します。

√
p2 + m1 2 + p
mK 2 − m1 2
p=
2mK

mK =

12

3.1

3 粒子衝突と崩壊 ∼ 静止 K + 崩壊を中心に ∼２体崩壊

例えば、Ke2 の場合単色運動量 pe = 246.8 MeV/c の e+ が、Kµ2 の場合 pµ = 235.5 MeV/c の µ+ である。一方。

Kπ2 の場合は m2 = 0 近似できないので、ちゃんと計算しましょう。
√√
mK = p2 + m1 2 + p2 + m2 2√
mK 2 − 2mK p2 + m2 2 + p2 + m2 2 = p2 + m1 2√
2mK p2 + m2 2 = m2K − m21 + m22( 2)2
mK − m21 + m22
2
2
p + m2 =2mK√()2mK 2 − m1 2 + m2 2p=
− m2 22mK
つまり、pπ = 205.4 MeV/c の π + が放出されるみたいな。m2 = 0 で前の式と一致してますね。

sample6.cxx

void sample6(){
double mK=493.677, me=0.511, mmu=105.6, mpi=139.6,mpi0=134.1;
double pK=0, pe=0,pnu=0, pmu=0, ppi=0, ppi0;
twobodydecay(mK, me, pe, 0, pnu);
twobodydecay(mK, mmu, pmu, 0, pnu);
twobodydecay(mK, mpi, ppi, mpi0, ppi0);
cout<<pe<<" "<<pmu<<" "<<ppi<<" "<<endl;
return;}void twobodydecay(double parentM,
double &daughterM1, double &daughterP1,
double &daughterM2, double &daughterP2){
daughterP1 = sqrt(((parentM**2 - daughterM1**2 + daughterM2**2)/(2*parentM))**2-daughterM2**2);
daughterP2 = daughterP1;
return;}

このプログラムは２体崩壊関数 twobodydecay() を使用して、娘粒子の運動量を出力します。引数に親粒子の質量、そして娘粒子の質量と運動量が返ってくる。sample6.cxx ソースの中身で最初の twobodydecay() は、Ke2 で e+ と

νe の２粒子の運動量を計算しています。引数の変数の前に & をつけると関数内で計算した値を、戻り値で返すことができるので便利ですよ。
親粒子が飛んでいる際意中に崩壊する事象を in-flight 崩壊と言うのはご存知かと。せっかくなのでローレンツ変換を使って、親粒子が in-flight 崩壊した時に放射される娘粒子の運動量ベクトルを決定していきましょう。π 0 は寿命
が 8.5 × 10−17 s なので、25 nm しか進まず、2 本のγ線を放出します。Kπ2 で生成された単色 π 0 が放出する 2γ を考えていきましょう。

sample7.cxx

void sample7(){double Eg1, pg1,Eg2, pg2,p[2][3]={0};double mpi=134.9, mg1=0,mg2=0;

double costheta,phi,sintheta,costhetagg;

double ppi[3]={0,0,205.4},mom;double beta[3],gamma[3];

13



3.1

3 粒子衝突と崩壊 ∼ 静止 K + 崩壊を中心に ∼２体崩壊

sample7.cxx の続きhist

gSystem->Load("sample6.cxx");

hist
2000

mom=0;
for(int i=0;i<3;i++) mom+=ppi[i];1800

Entries
Mean
RMS

10000
0.12540.3315

0.8

1

1600

for(int i=0;i<3;i++){1400

beta[i]=pK[i]/sqrt(mK**2 + momK**2);

1200

gamma[i]=1./sqrt(1.-beta[i]**2);

1000
800

}

600
400

hist=new TH1F("hist","hist",100,-1,1);

200
0
−1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

twobodydecay(mpi, mg1, pg1, mg2, pg2);
Eg1=sqrt(mg1**2+pg1**2);
Eg2=sqrt(mg2**2+pg2**2);;
for(int i=0;i<10000;i++){costheta=gRandom->Uniform(-1,1);
phi=2*3.141592*gRandom->Uniform(0,1);
sintheta=sqrt(1-costheta**2);
p[0][0]=pg1*sintheta*sin(phi);p[0][1]=pg1*sintheta*cos(phi);p[0][2]=pg1*costheta;
p[1][0]=-pg2*sintheta*sin(phi);p[1][1]=-pg2*sintheta*cos(phi);p[1][2]=-pg2*costheta;for(int j=0;j<3;j++){p[0][j] = gamma[j]*beta[j]*Emu + gamma[j]*p[0][j];p[1][j] = gamma[j]*beta[j]*Enu + gamma[j]*p[1][j];}costhetagg = (p[0][0]*p[1][0]+p[0][1]*p[1][1]+p[0][2]*p[1][2])/sqrt(p[1][0]**2+p[1][1]**2+p[1][2]**2)/sqrt(p[0][0]**2+p[0][1]**2+p[0][2]**2);
hist->Fill(costhetagg);}
c1=new TCanvas("c1");
hist->Draw();
c1->Print("output.pdf");
return;
}




14

3.2

3 粒子衝突と崩壊 ∼ 静止 K + 崩壊を中心に ∼3 体崩壊

このことからも、π 0 から放出されたγ線の opening 角度はローレンツブーストによって前方方向に伸ばされていることがわかります。gSystem->Load() を実行すると、sample6.cxx の関数を読み込むことができます。cos θ = 0.4
は角度でいうと 66.4 度くらい。-1 までテイルをひくのは π 0 静止系で等方的にγ線が放射する中で丁度 π 0 の進
行方向と平行な事象が存在するから。0.4 にピークを持つのは、γ線の向きと π 0 の向きが垂直の時の角度相関

cos(2 tan−1 (134.9/205.4)) = 0.397 が最大であるからで、等方的に放射してもこれ以上角度が狭くなることが理論的にないので、絶壁を持った連続分布を生成します。

3.2 3 体崩壊
主な 3 体崩壊は Ke3 , Kµ3 , K + → π + π + π − , K + → π + π 0 π 0 などの π を伴います。また、厄介なのはそれぞれの娘粒子の運動量は連続的に分布することです。エネルギー・運動量保存則を満たすように幾何学的な要請とその中での崩壊密度分布から、娘粒子の運動量が決定されます。この分布を Dalitz plot といいます。

3.2.1 K → 3π
まずは K → 3π 過程から見ていきましょう。π− 静止系では π + の 2 体崩壊と同等。その時の π + の運動量は√mK = 2 mπ 2 + p 2 + mπ√
mK − mπ = 2 mπ 2 + p2
mK − mπ= mπ 2 + p2
2√()2
mK − mπ− mπ 2p=2(2)で表せる。3 体崩壊は連続的な運動量を持ち、π + の最大の運動量は 110.15 MeV であると計算できました。K → 3πの Dalitz plot 分布は以下の散乱振幅から表されます。

{}2s3 − s0s3 − s0|M | ∝ 1 + g
+h
m2π+
m2π+
{
}2
s2 − s1
s2 − s1
+
k
+ ...
+j
m2π+
m2π+
2

(3)

ここで、mπ+ は π + の質量、

si = (PK − Pi )2= (mK − mi )2 − 2mK Ti , (i = 1, 2, 3),
1∑
s0 =
si
3 i
=

1
(mK 2 + m1 2 + m2 2 + m3 2 )
3

(4)

(5)

ここで、Pi と mi , Ti はそれぞれ崩壊後の娘核 i の 4 元運動量と質量、運動エネルギーを示す。パラメータ g, h, k
は最新の研究結果を採用している。CP 対称性が破れていないと仮定しているため j = 0. 粒子 3 は奇数パイオンを示しています。K ± → π ± π ± π ∓ において

g = −0.21134 ± 0.00017,−2

h = (1.848 ± 0.040) × 10
j = 0,k = (−4.63 ± 0.14) × 10−3です。この時の Dalitz plot を描いて見ましょう。

15

(6)
,

(7)
(8)
(9)

3.2

3 粒子衝突と崩壊 ∼ 静止 K + 崩壊を中心に ∼

3 体崩壊


sample8.cxx

void sample8(){gStyle->SetOptStat(0);double mK=493.677, mpi=139;
double p[3][3]={0};dalitz23=new TH2F("dalitz23","dalitz23",120,0,120,120,0,120);
dalitz31=new TH2F("dalitz31","dalitz31",120,0,120,120,0,120);
dalitz12=new TH2F("dalitz12","dalitz12",120,0,120,120,0,120);for(int i=0;i<100000;i++){
DalitzdecayKto3piPlus(mK,mpi,p[0][0],p[0][1],p[0][2],mpi,p[1][0],p[1][1],p[1][2],mpi,p[2][0],p[2][1],p[2][2]);dalitz23->Fill(sqrt(p[1][0]**2+p[1][1]**2+p[1][2]**2),
sqrt(p[2][0]**2+p[2][1]**2+p[2][2]**2));
dalitz31->Fill(sqrt(p[2][0]**2+p[2][1]**2+p[2][2]**2),
sqrt(p[0][0]**2+p[0][1]**2+p[0][2]**2));
dalitz12->Fill(sqrt(p[0][0]**2+p[0][1]**2+p[0][2]**2),
sqrt(p[1][0]**2+p[1][1]**2+p[1][2]**2));}
c1=new TCanvas("c1","",800,300);
c1->Divide(3,1);
c1->cd(1);dalitz23->SetXTitle("p_{#pi2} (MeV/c)");dalitz23->GetXaxis()->CenterTitle();dalitz23->SetYTitle("p_{#pi3} (MeV/c)");dalitz23->GetYaxis()->CenterTitle();dalitz23->Draw("colz");c1->cd(2);dalitz31->SetXTitle("p_{#pi3} (MeV/c)");dalitz31->GetXaxis()->CenterTitle();dalitz31->SetYTitle("p_{#pi1} (MeV/c)");dalitz31->GetYaxis()->CenterTitle();dalitz31->Draw("colz");c1->cd(3);dalitz12->SetXTitle("p_{#pi1} (MeV/c)");dalitz12->GetXaxis()->CenterTitle();dalitz12->SetYTitle("p_{#pi2} (MeV/c)");dalitz12->GetYaxis()->CenterTitle();dalitz12->Draw("colz");
c1->Print("output.pdf");
return;
}



16

3.2

3 粒子衝突と崩壊 ∼ 静止 K + 崩壊を中心に ∼

3 体崩壊


sample8.cxx の続きvoid DalitzdecayKto3piPlus(
double parentM,
double daughter1M,

double &daughter1P1,double &daughter1P2,

double &daughter1P3,double daughter2M,

double &daughter2P1,

double &daughter2P2,double &daughter2P3,double daughter3M,
double &daughter3P1,double &daughter3P2,double &daughter3P3){
double en_p[3]={0},en_E[3]={0},s[4]={0},en_r;
double pmax=sqrt(((parentM - daughter1M)/2)**2 - daughter1M**2);
double Tmax=(parentM - 3*daughter1M)/2;
double Emax=(parentM - daughter1M)/2;
double g=-0.21134;
double h=1.848e-2;
double k=-4.63e-3;
double M, M2;
double costheta, sintheta, phi, sintheta3decay, phin, costheta3decay;do{en_E[0] = gRandom->Uniform(daughter1M,Tmax+daughter1M);
en_E[1] = gRandom->Uniform(daughter1M,Tmax+daughter1M);
en_E[2] = parentM - en_E[0] - en_E[1];for(int i=0;i<3;i++) en_p[i]=sqrt(en_E[i]**2 - daughter1M**2);s[0]=(parentM**2 + 3*daughter1M**2)/3;
for(int i=1;i<4;i++){s[i]=(parentM-daughter1M)**22*parentM*(en_E[i-1]-daughter1M);}
M = en_E[0]+en_E[1]+en_E[2];M2 = 1. + g*(s[3]-s[0])/daughter1M**2+ h*(s[3]-s[0])**2/(daughter1M**4)+ k*(s[2]-s[1])**2/(daughter1M**4);
M2*=1e-1;}while(en_E[0] > Tmax+daughter1M ||en_E[1] > Tmax+daughter1M ||en_E[2] > Tmax+daughter1M ||en_E[0] < daughter1M ||en_E[1] < daughter1M ||en_E[2] < daughter1M ||gRandom->Uniform(0.00,1.00) > M2
);




17

3.2

3 粒子衝突と崩壊 ∼ 静止 K + 崩壊を中心に ∼

3 体崩壊


sample8.cxx の続きcostheta3decay =(en_E[2]**2 - en_E[0]**2 - en_E[1]**2 + daughter1M**2)/(2*en_p[0]*en_p[1]);// pion 1
costheta = -1. +2.*gRandom->Uniform(0,1);
sintheta = sqrt((1.-costheta)*(1.+costheta));phi = 2.*3.141592*gRandom->Uniform(0,1);
daughter1P1 = en_p[0]*sintheta*cos(phi);daughter1P2 = en_p[0]*sintheta*sin(phi);daughter1P3 = en_p[0]*costheta;// pion 2
sintheta3decay=sqrt((1.+costheta3decay)*(1.-costheta3decay));phin=2.*3.141592*gRandom->Uniform(0,1);
daughter2P1 = en_p[1]*(sintheta3decay*cos(phin)*costheta*cos(phi)- sintheta3decay*sin(phin)*sin(phi)+ costheta3decay*sintheta*cos(phi));daughter2P2 = en_p[1]*(sintheta3decay*cos(phin)*costheta*sin(phi)+ sintheta3decay*sin(phin)*cos(phi)+ costheta3decay*sintheta*sin(phi));daughter2P3 = en_p[1]*(- sintheta3decay*cos(phin)*sintheta+ costheta3decay*costheta);// pion 3
daughter3P1 = - daughter1P1 - daughter2P1;
daughter3P2 = - daughter1P2 - daughter2P2;
daughter3P3 = - daughter1P3 - daughter2P3;
return;
}




18

3.2

3 粒子衝突と崩壊 ∼ 静止 K + 崩壊を中心に ∼3 体崩壊

DalitzdecayKto3piPlus() 関数の中では運動学的に３つの π のエネルギーを与えて、散乱振幅式 (3) から崩壊強度を決定しています。２粒子の運動量がわかると、その間の角度を計算できる。

p1 · p2|p1 | |p2 |(p1 + p2 )2 − p1 2 − p2 2=
2 |p1 | |p2 |cos Θ ==

(p3 )2 − (E1 2 − m1 2 ) − (E2 2 − m2 2 )2 |p1 | |p2 |=(E3 2 − m3 2 ) − (E1 2 − m1 2 ) − (E2 2 − m2 2 )2 |p1 | |p2 |(10)

ここで m1 = m2 = m3 = mπ なのでキャンセルして、

cos Θ =

E3 2 − E1 2 − E2 2 + mπ 2
2 |p1 | |p2 |(11)と計算できます。次に３粒子の運動量の成分を決定しましょう。粒子１は θ, ϕ の一様乱数を使って
sin θ cos ϕ
p1 = |p1 | ·  sin θ sin ϕ 
cos θ

(12)

とわかる。粒子２は粒子 1 の周りを角度 Θ で回転する自由度を持つ。一様乱数 Φ を用いて、


sin Θ cos Φ cos θ cos ϕ − sin Θ sin Φ sin ϕ + cos Θ sin Θ cos ϕ
p2 = |p2 | ·  sin Θ cos Φ cos θ sin ϕ + sin Θ sin Φ cos ϕ + cos Θ sin Θ sin ϕ − sin Θ cos Φ sin θ + cos Θ cos θ(13)と計算することができる。粒子 3 は粒子 1, 2 と親核の運動量保存則から

p3 = −p1 − p2(14)と計算することができる。
dalitz23

dalitz31

dalitz12

120
70

80

50

80

60

50

p π 1 (MeV/c)100p π 3 (MeV/c)

60

60

40

60

30

40

30

40
20

20

0
0

10

20

40

図2

60
80
p π 2 (MeV/c)

100

120

0

60

100

50

80

40

70

100

40

120

p π 2 (MeV/c)

120

60

30

40

20

20

10

20
20

0
0

10

20

40

60
80
p π 3 (MeV/c)

100

120

0

0
0

20

40

60
80
p π 1 (MeV/c)

100

120

0

K + → π + π + π − における π の運動量分布. 粒子 1, 2 は π + , 粒子 3 は π − を示す.

K ± → π ± π 0 π 0 において
g0 = 0.626 ± 0.007,
h0 = 0.052 ± 0.008,

(15)
(16)

j0 = 0,k0 = 0.0054 ± 0.0035

(17)
(18)

のパラメータで同様に計算できる。

19

3.2

3 粒子衝突と崩壊 ∼ 静止 K + 崩壊を中心に ∼

3 体崩壊


Dalitzdecayto3piZero.cxx

void DalitzdecayKto3piZero(double parentM,
double daughter1M,

double &daughter1P1,double &daughter1P2,

double &daughter1P3,double daughter2M,

double &daughter2P1,

double &daughter2P2,double &daughter2P3,double daughter3M,
double &daughter3P1,double &daughter3P2,double &daughter3P3){
double en_p[3]={0},en_E[3]={0},s[4]={0},en_r;
double Tmax=(parentM - daughter1M - daughter2M - daughter3M)/2;
double g=0.626;
double h=0.052;
double k=0.0054;
double M, M2;
double costheta, sintheta, phi, sintheta3decay, phin, costheta3decay;do{en_E[0] = gRandom->Uniform(daughter1M,Tmax+daughter1M);
en_E[1] = gRandom->Uniform(daughter2M,Tmax+daughter2M);
en_E[2] = parentM - en_E[0] - en_E[1];en_p[0]=sqrt(en_E[0]**2 - daughter1M**2);
en_p[1]=sqrt(en_E[1]**2 - daughter2M**2);
en_p[2]=sqrt(en_E[2]**2 - daughter3M**2);s[0]=(parentM**2 + daughter1M**2+daughter2M**2+daughter3M**2)/3;s[1]=(parentM-daughter1M)**2-2*parentM*(en_E[0]-daughter1M);
s[2]=(parentM-daughter2M)**2-2*parentM*(en_E[1]-daughter2M);
s[3]=(parentM-daughter3M)**2-2*parentM*(en_E[2]-daughter3M);
M = en_E[0]+en_E[1]+en_E[2];M2 = 1. + g*(s[3]-s[0])/daughter3M**2+ h*(s[3]-s[0])**2/(daughter3M**4)+ k*(s[2]-s[1])**2/(daughter3M**4);
M2*=1e-1;}while(en_E[0] > Tmax+daughter1M ||en_E[1] > Tmax+daughter2M ||en_E[2] > Tmax+daughter3M ||en_E[0] < daughter1M ||en_E[1] < daughter2M ||en_E[2] < daughter3M ||gRandom->Uniform(0.00,1.00) > M2
);



20

3.2

3 粒子衝突と崩壊 ∼ 静止 K + 崩壊を中心に ∼

3 体崩壊


Dalitzdecayto3piZero.cxx の続きcostheta3decay =(en_E[2]**2 - daughter3M**2 - en_E[0]**2 + daughter1M**2 - en_E[1]**2+ daughter2M**2)/(2*en_p[0]*en_p[1]);costheta = -1. +2.*gRandom->Uniform(0,1);// pion 1

sintheta = sqrt((1.-costheta)*(1.+costheta));phi = 2.*3.141592*gRandom->Uniform(0,1);
daughter1P1 = en_p[0]*sintheta*cos(phi);daughter1P2 = en_p[0]*sintheta*sin(phi);daughter1P3 = en_p[0]*costheta;
sintheta3decay=sqrt((1.+costheta3decay)*(1.-costheta3decay));// pion 2

phin=2.*3.141592*gRandom->Uniform(0,1);
daughter2P1 = en_p[1]*(sintheta3decay*cos(phin)*costheta*cos(phi)- sintheta3decay*sin(phin)*sin(phi)+ costheta3decay*sintheta*cos(phi));daughter2P2 = en_p[1]*(sintheta3decay*cos(phin)*costheta*sin(phi)+ sintheta3decay*sin(phin)*cos(phi)+ costheta3decay*sintheta*sin(phi));daughter2P3 = en_p[1]*(- sintheta3decay*cos(phin)*sintheta+ costheta3decay*costheta);
daughter3P1 = - daughter1P1 - daughter2P1;// pion 3

daughter3P2 = - daughter1P2 - daughter2P2;
daughter3P3 = - daughter1P3 - daughter2P3;
return;}dalitz23

dalitz31

dalitz12

120

120

60
100

100

40

60

p π 1 (MeV/c)

p π 3 (MeV/c)

80
40

30

60

30

40

90
100

50
80

100

50

80
70

80

60

p π 2 (MeV/c)

120

60

50
40

40

20

20

10

20

0
0

0

0
0

40

30

20
20

10

20
10

0
0

20

40

60
80
p π 2 (MeV/c)

100

120

0

20

40

60
80
p π 3 (MeV/c)

100

120

20

40

60
80
p π 1 (MeV/c)100

図 3 K + → π 0 π 0 π + における π の運動量分布. 粒子 1, 2 は π 0 , 粒子 3 は π + を示す.

21

120

0

3.2

3 粒子衝突と崩壊 ∼ 静止 K + 崩壊を中心に ∼3 体崩壊

3.2.2 K → πlνν 静止系において K → πlν の娘粒子である l と π 0 の運動量は等しく、この運動量が最大値をとる。mν = 0 を仮定して、Ke3 の場合、

mK =m2K − 2mk√√me 2 + p2 + mπ0 2 + p2√
mπ0 2 + p2 + mπ0 2 + p2 = me 2 + p2
m2 − mπ0 2 + me 2
mπ0 2 + p2 = k2mk√2
mk − mπ0 2 + me 2p=
− mπ 0 22mk∼ 228.7 MeV/c.

(19)Kµ3 の場合、p ∼ 215.5 MeV/c.

(20)

３体崩壊を解くにおいて運動量の取りうる範囲をまず制限します。その次に確率密度から、崩壊頻度を決定します。

K → πlν (Ke3 , Kµ3 ) の確率密度は以下で表せます。
{}
2
ρ(t) ∝ f (t) A + B ξ(t) + C ξ (t)2(21)

ここで、使用されているパラメータはt = mK 2 + mπ0 2 − 2mK Eπ ,()tξ(t) = ξ0 1 + λ,mπ 0 2tf (t) = 1 + λ,mπ 0 2
mK 2 + m2π0 − mlE=
− Eπ 02mKA = mK (2El Eν − mK E) + ml 2()
E2
B = ml Eν −2
ml 2 EC=
4です。
ここで f (t) を Form Factor (形状因子) という。

Ke3 の場合、初期値として λ = 0.0286, ξ0 = −0.35.
Kµ3 の場合、λ = 0.033, ξ0 = −0.35.

22

(22)
(23)
(24)
(25)
(

E
− Eν
4

)
(26)
(27)
(28)

3.2

3 粒子衝突と崩壊 ∼ 静止 K + 崩壊を中心に ∼

3 体崩壊


sample9.cxx

void sample9(){gStyle->SetOptStat(0);double mK=493.677, mmu=105.6, mpi=139,mpi0=134,me=0.511,mnu=0;
double pK[3]={0};dalitz23=new TH2F("dalitz23","dalitz23",100,0,250,100,0,300);
dalitz31=new TH2F("dalitz31","dalitz31",100,0,250,120,0,300);
for(int i=0;i<100000;i++){
DalitzdecayKl3(mK,mpi0,p[0][0],p[0][1],p[0][2],me,p[1][0],p[1][1],p[1][2],mnu,p[2][0],p[2][1],p[2][2]);dalitz23->Fill(sqrt(me**2 + p[1][0]**2+p[1][1]**2+p[1][2]**2),
sqrt(mpi0**2+ p[0][0]**2+p[0][1]**2+p[0][2]**2));
DalitzdecayKl3(mK,mpi0,p[0][0],p[0][1],p[0][2],mmu,p[1][0],p[1][1],p[1][2],mnu,p[2][0],p[2][1],p[2][2]);dalitz31->Fill(sqrt(mmu**2+ p[1][0]**2+p[1][1]**2+p[1][2]**2),
sqrt(mpi0**2+ p[0][0]**2+p[0][1]**2+p[0][2]**2));if(i%100==0)cout<<"."<<flush;}
c1=new TCanvas("c1","",700,300);
c1->Divide(2,1);
c1->cd(1);dalitz23->SetTitle("K^{+}_{e3}");
dalitz23->SetXTitle("Electron Energy (MeV)");
dalitz23->GetXaxis()->CenterTitle();dalitz23->SetYTitle("Pion Energy (MeV)");
dalitz23->GetYaxis()->CenterTitle();dalitz23->Draw("colz");




23

3.2

3 粒子衝突と崩壊 ∼ 静止 K + 崩壊を中心に ∼

3 体崩壊


sample9.cxx の続き

c1->cd(2);dalitz31->SetTitle("K^{+}_{#mu3}");
dalitz31->SetXTitle("Muon Energy (MeV)");
dalitz31->GetXaxis()->CenterTitle();dalitz31->SetYTitle("Pion Energy (MeV)");
dalitz31->GetYaxis()->CenterTitle();dalitz31->Draw("colz");
c1->Print("output.pdf");
return;}void DalitzdecayKl3(
double parentM,
double daughter1M,double &daughter1P1,double &daughter1P2,double &daughter1P3,
double daughter2M,double &daughter2P1,double &daughter2P2,double &daughter2P3,
double daughter3M,double &daughter3P1,double &daughter3P2,double &daughter3P3){
double en_p[3]={0},en_E[3]={0};double pmax=
sqrt(((parentM**2 + daughter1M**2 daughter2M**2)/2/parentM)**2 - daughter1M**2);
double M;
double costheta, sintheta, phi, sintheta3decay, phin, costheta3decay;
double Epi,Epi_max,E,q2,F,Fmax,Xi,coeffA,coeffB,coeffC,RhoMax,Rho,Epi,El,Enu;double massK= parentM;

double massPi = daughter1M;
double massL= daughter2M;

double massNu = daughter3M;
double pLambda, pXi0;
if(massL==0.511){ pLambda = 0.0286; pXi0 = -0.35;}//0.0297,-0.3
else if (massL=105.6){pLambda = 0.033;pXi0 = -0.35;}//0.0297,0.17
do{
en_E[0] = gRandom->Uniform(daughter1M,sqrt(daughter1M**2+pmax**2));
en_E[1] = gRandom->Uniform(daughter2M,sqrt(daughter2M**2+pmax**2));
en_E[2] = parentM - en_E[0] - en_E[1];//nu
en_p[0]=sqrt(en_E[0]**2 - daughter1M**2);
en_p[1]=sqrt(en_E[1]**2 - daughter2M**2);
en_p[2]=sqrt(en_E[2]**2 - daughter3M**2);



24

3.2

3 粒子衝突と崩壊 ∼ 静止 K + 崩壊を中心に ∼

3 体崩壊


sample9.cxx の続き

Epi=en_E[0];El =en_E[1];Enu=en_E[2];Epi_max = (massK*massK+massPi*massPi-massL*massL)/2.0/massK;E= Epi_max - Epi;

M = en_E[0]+en_E[1]+en_E[2];q2 = massK*massK + massPi*massPi - 2.0*massK*Epi;
F = 1.0 + pLambda*q2/massPi/massPi;
Fmax = 1.0;
if (pLambda >0.0) Fmax = (1.0 + pLambda*(massK*massK/massPi/massPi+1.0));
Xi = pXi0*(1.0 + pLambda*q2/massPi/massPi);
coeffA = massK*(2.0*El*Enu-massK*E)+massL*massL*(E/4.0-Enu);
coeffB = massL*massL*(Enu-E/2.0);
coeffC = massL*massL*E/4.0;
RhoMax = (Fmax*Fmax)*(massK*massK*massK/8.0);
Rho = (F*F)*(coeffA + coeffB*Xi + coeffC*Xi*Xi);
M=en_E[0]+en_E[1]+en_E[2];costheta3decay =(en_E[2]**2 - daughter3M**2 - en_E[0]**2+ daughter1M**2 - en_E[1]**2 + daughter2M**2)/(2*en_p[0]*en_p[1]);}while(en_E[0] > sqrt(daughter1M**2+pmax**2) ||en_E[1] > sqrt(daughter2M**2+pmax**2) ||en_E[2] > sqrt(daughter3M**2+pmax**2) ||en_E[0] < daughter1M ||en_E[1] < daughter2M ||en_E[2] < daughter3M ||costheta3decay>1 ||costheta3decay<-1||
gRandom->Uniform(0.00,1.00) > Rho/RhoMax||
Rho/RhoMax < 0
);




25

3.2

3 粒子衝突と崩壊 ∼ 静止 K + 崩壊を中心に ∼

3 体崩壊


sample9.cxx の続きcostheta = -1. +2.*gRandom->Uniform(0,1);// particle 1

sintheta = sqrt((1.-costheta)*(1.+costheta));phi = 2.*3.141592*gRandom->Uniform(0,1);
daughter1P1 = en_p[0]*sintheta*cos(phi);daughter1P2 = en_p[0]*sintheta*sin(phi);daughter1P3 = en_p[0]*costheta;

sintheta3decay=sqrt((1.+costheta3decay)*(1.-costheta3decay));// particle 2

phin=2.*3.141592*gRandom->Uniform(0,1);
daughter2P1 = en_p[1]*(sintheta3decay*cos(phin)*costheta*cos(phi)- sintheta3decay*sin(phin)*sin(phi)+ costheta3decay*sintheta*cos(phi));daughter2P2 = en_p[1]*(sintheta3decay*cos(phin)*costheta*sin(phi)+ sintheta3decay*sin(phin)*cos(phi)+ costheta3decay*sintheta*sin(phi));daughter2P3 = en_p[1]*(- sintheta3decay*cos(phin)*sintheta+ costheta3decay*costheta);
daughter3P1 = - daughter1P1 - daughter2P1;// particle 3

daughter3P2 = - daughter1P2 - daughter2P2;
daughter3P3 = - daughter1P3 - daughter2P3;
return;
}



+

+

Kµ3

K e3
260
260

120

100
240

240

100
220

Pion Energy (MeV)

Pion Energy (MeV)

80

220

80

200

200

60

60

180

180
40

160

40

160
140

20

140

20
120

120
0

50

100
150
Electron Energy (MeV)

図4

200

250

0

100

120

140
160
180
200
Muon Energy (MeV)

Kl3 の Dalitz Plot. Ke3 (左), Kµ3 (右)

26

220

240

0

3.3

3 粒子衝突と崩壊 ∼ 静止 K + 崩壊を中心に ∼静止 K 崩壊荷電粒子のスペクトラム

3.3 静止 K 崩壊荷電粒子のスペクトラム
荷電 K + 中間子が崩壊した時のほとんどの崩壊過程をみて来ました。磁場を用いた実験で静止 K + から崩壊した荷電粒子は (1) プラスチックシンチレータ検出器でできた TOF start を通過し、(2) ドリフトチェンバーを通過し、磁場で曲がって (3) ドリフトチェンバーを通過し、(4) プラスチックシンチレータ検出器でできた TOF stop を通過するセットアップを考えます。ここでは、簡単のためトラッキングによる運動量分解能が 0.9%、TOF の時間分解能
が 100 ps であると仮定して、磁場での軌道は 10 m に規格化されているとする。

Ke2 と Ke3 では e+ が、Kµ2 と Kµ3 では µ+ が、そして Kπ2 と K + → π + π + π − , K + → π 0 π 0 π + では π + 放出
される。これらの粒子の得られる運動量を描いてみましょう。

sample10.cxx

void sample10(){gStyle->SetOptStat(0);double mK=493.677, mmu=105.6, mpi=139,mpi0=134,me=0.511,mnu=0;
double pmu,pnu,pe,ppi,ppi0,p[3][3]={0};double costheta,phi,sintheta,rndm;
double Gamma1=1.58e-5; //Ke2
double Gamma2=63.55e-2;//Kmu2
double Gamma3=5.07e-2; //Ke3
double Gamma4=3.35e-2; //Kmu3
double Gamma5=20.66e-2;//Kpi2
double Gamma6=1.76e-2; //Kpi3plus
double Gamma7=5.59e-2; //Kpi3zero
double dmy;
double resoP=0.9e-2;
total=new TH1F("total","total",600,0,300);
Ke2=new TH1F("Ke2","Ke2",600,0,300);
Kmu2=new TH1F("Kmu2","Kmu2",600,0,300);
Ke3=new TH1F("Ke3","Ke3",600,0,300);
Kmu3=new TH1F("Kmu3","Kmu3",600,0,300);
Kpi2=new TH1F("Kpi2","Kpi2",600,0,300);
Kpi3plus=new TH1F("Kpi3plus","Kpi3plus",600,0,300);
Kpi3zero=new TH1F("Kpi3zero","Kpi3zero",600,0,300);
gSystem->Load("sample6.cxx");
gSystem->Load("sample8.cxx");
gSystem->Load("DalitzdecayKto3piZero.cxx");
gSystem->Load("sample9.cxx");




27

3.3

3 粒子衝突と崩壊 ∼ 静止 K + 崩壊を中心に ∼静止 K 崩壊荷電粒子のスペクトラム

sample10.cxx の続きfor(int i=0;i<10000000;i++){rndm=gRandom->Uniform(0,1);
if(rndm<=Gamma1){
twobodydecay(mK, me, pe, 0, pnu);
dmy=gRandom->Gaus(pe,pe*resoP);
total->Fill(dmy);
Ke2->Fill(dmy);}else if(rndm>Gamma1 && rndm<=Gamma1+Gamma2){
twobodydecay(mK, mmu, pmu, 0, pnu);
dmy=gRandom->Gaus(pmu,pmu*resoP);
total->Fill(dmy);
Kmu2->Fill(dmy);}else if(rndm>Gamma1+Gamma2&& rndm<=Gamma1+Gamma2+Gamma3){
DalitzdecayKl3(mK,mpi0,p[0][0],p[0][1],p[0][2],me,p[1][0],p[1][1],p[1][2],mnu,p[2][0],p[2][1],p[2][2]);pe=sqrt(p[1][0]**2+p[1][1]**2+p[1][2]**2);
dmy=gRandom->Gaus(pe,pe*resoP);
total->Fill(dmy);
Ke3->Fill(dmy);}else if(rndm>Gamma1+Gamma2+Gamma3&& rndm<=Gamma1+Gamma2+Gamma3+Gamma4){
DalitzdecayKl3(mK,mpi0,p[0][0],p[0][1],p[0][2],mmu,p[1][0],p[1][1],p[1][2],mnu,p[2][0],p[2][1],p[2][2]);pmu=sqrt(p[1][0]**2+p[1][1]**2+p[1][2]**2);
dmy=gRandom->Gaus(pmu,pmu*resoP);
total->Fill(dmy);
Kmu3->Fill(dmy);}else if(rndm>Gamma1+Gamma2+Gamma3+Gamma4&&rndm<=Gamma1+Gamma2+Gamma3+Gamma4+Gamma5){
twobodydecay(mK, mpi, ppi, mpi0, ppi0);
dmy=gRandom->Gaus(ppi,ppi*resoP);
total->Fill(dmy);
Kpi2->Fill(dmy);


}


28

3.3

3 粒子衝突と崩壊 ∼ 静止 K + 崩壊を中心に ∼静止 K 崩壊荷電粒子のスペクトラム

sample10.cxx の続きelse if(rndm>Gamma1+Gamma2+Gamma3+Gamma4+Gamma5&& rndm<=Gamma1+Gamma2+Gamma3+Gamma4+Gamma5+Gamma6){
DalitzdecayKto3piPlus(mK,mpi,p[0][0],p[0][1],p[0][2],mpi,p[1][0],p[1][1],p[1][2],mpi,p[2][0],p[2][1],p[2][2]);ppi=sqrt(p[0][0]**2+p[0][1]**2+p[0][2]**2);
dmy=gRandom->Gaus(ppi,ppi*resoP);
total->Fill(dmy);
Kpi3plus->Fill(dmy);}else if(rndm>Gamma1+Gamma2+Gamma3+Gamma4+Gamma5+Gamma6&& rndm<=Gamma1+Gamma2+Gamma3+Gamma4+Gamma5+Gamma6+Gamma7){
DalitzdecayKto3piZero(mK,
mpi0,p[0][0],p[0][1],p[0][2],mpi0,p[1][0],p[1][1],p[1][2],mpi,p[2][0],p[2][1],p[2][2]);ppi=sqrt(p[2][0]**2+p[2][1]**2+p[2][2]**2);
dmy=gRandom->Gaus(ppi,ppi*resoP);
total->Fill(dmy);
Kpi3zero->Fill(dmy);}
if(i%10000==0)cout<<"."<<flush;}
c1=new TCanvas("c1","",500,300);
c1->SetLogy();
total->SetTitle("");
total->SetXTitle("Momenum p (MeV/c)");total->GetXaxis()->CenterTitle();total->SetLineColor(1);total->SetFillColor(1);total->SetFillStyle(3003);
total->Draw();
Ke2->SetLineColor(2);Ke2->Draw("same");Kmu2->SetLineColor(3);Kmu2->Draw("same");Ke3->SetLineColor(4);Ke3->Draw("same");Kmu3->SetLineColor(5);Kmu3->Draw("same");



29

3.3

3 粒子衝突と崩壊 ∼ 静止 K + 崩壊を中心に ∼静止 K 崩壊荷電粒子のスペクトラム

sample10.cxx の続き

Kpi2->SetLineColor(6);Kpi2->Draw("same");
Kpi3plus->SetLineColor(7);Kpi3plus->Draw("same");
Kpi3zero->SetLineColor(8);Kpi3zero->Draw("same");c1->Print("output.pdf");
return;}

２体崩壊における荷電粒子のスペクトラムはピーク構造なので、Ke2 , Kµ2 , Kπ2 はわかりやすいですね。3 体崩壊で放出される荷電粒子は連続的な運動量を持つことが視覚的にも理解できると思います。Kµ2 の裾に Ke2 が隠れているため、運動量分布から Ke2 を観測することは難しいことがわかります。e/µ の粒子識別によって Ke2 ピークを見やすくできるはずです。

6

10

5

10

4

10

3

10

2

10

10

1
0

図5

50

100

150
Momenum p (MeV/c)

200

250

300

静止 K + 崩壊 運動量スペクトラム. Ke2 (赤), Kµ2 (緑), Ke3 (青), Kµ3 (黄), Kπ2 (マゼンタ), K + →π + π + π − (暗緑), K + → π 0 π 0 π + (シアン).

TOF を用いて e+ /µ+ /pi+ の識別をしていきます。運動量 p =200–250 MeV/c の速度比 β はそれぞれ
βe = 1 − (2.089 ∼ 3.264) × 10−6βµ = 0.884 ∼ 0.921
βπ = 0.821 ∼ 0.874なので、10 m 走ると t = L/β だけ時間がかかる。このことから、運動量と時間差の分布を解析していきましょう。

30

3.3

3 粒子衝突と崩壊 ∼ 静止 K + 崩壊を中心に ∼静止 K 崩壊荷電粒子のスペクトラム

sample11.cxx

void sample11(){
gStyle->SetOptStat(0);double mK=493.677, mmu=105.6, mpi=139,mpi0=134,me=0.511,mnu=0;
double pmu,pnu,pe,ppi,ppi0,p[3][3]={0},P;
double costheta,phi,sintheta,rndm;
double beta,c=3.0e8,dt;
double Gamma1=1.58e-5; //Ke2
double Gamma2=63.55e-2;//Kmu2
double Gamma3=5.07e-2; //Ke3
double Gamma4=3.35e-2; //Kmu3
double Gamma5=20.66e-2;//Kpi2
double Gamma6=1.76e-2; //Kpi3plus
double Gamma7=5.59e-2; //Kpi3zero
double dmy,dmy2;
double resoP=0.9e-2;//%
double resoT=0.1;//ns
file=new TFile("kaon.root","recreate");total=new TH2F("total","total",300,0,300,300,20,120);
kaon=new TTree("kaon","kaon");kaon->Branch("P",&P,"P/D");
kaon->Branch("dt",&dt,"dt/D");
gSystem->Load("sample6.cxx");
gSystem->Load("sample8.cxx");
gSystem->Load("DalitzdecayKto3piZero.cxx");
gSystem->Load("sample9.cxx");for(int i=0;i<10000000;i++){rndm=gRandom->Uniform(0,1);
if(rndm<=Gamma1){
twobodydecay(mK, me, pe, 0, pnu);
P=gRandom->Gaus(pe,pe*resoP);
dmy2=10/c/(pe/sqrt(me**2 +pe**2))*1e9;//ns
dt=gRandom->Gaus(dmy2,dmy2*resoT);
total->Fill(P,dt);
}




31

3.3

3 粒子衝突と崩壊 ∼ 静止 K + 崩壊を中心に ∼静止 K 崩壊荷電粒子のスペクトラム

sample11.cxx の続きelse if(rndm>Gamma1 && rndm<=Gamma1+Gamma2){
twobodydecay(mK, mmu, pmu, 0, pnu);
P=gRandom->Gaus(pmu,pmu*resoP);
dmy2=10/c/(pmu/sqrt(mmu**2 +pmu**2))*1e9;//ns
dt=gRandom->Gaus(dmy2,dmy2*resoT);
total->Fill(P,dt);}
else if(rndm>Gamma1+Gamma2&& rndm<=Gamma1+Gamma2+Gamma3){
DalitzdecayKl3(mK,mpi0,p[0][0],p[0][1],p[0][2],me,p[1][0],p[1][1],p[1][2],mnu,p[2][0],p[2][1],p[2][2]);pe=sqrt(p[1][0]**2+p[1][1]**2+p[1][2]**2);
P=gRandom->Gaus(pe,pe*resoP);
dmy2=10/c/(pe/sqrt(me**2 +pe**2))*1e9;//ns
dt=gRandom->Gaus(dmy2,dmy2*resoT);
total->Fill(P,dt);}
else if(rndm>Gamma1+Gamma2+Gamma3&& rndm<=Gamma1+Gamma2+Gamma3+Gamma4){
DalitzdecayKl3(mK,mpi0,p[0][0],p[0][1],p[0][2],mmu,p[1][0],p[1][1],p[1][2],mnu,p[2][0],p[2][1],p[2][2]);pmu=sqrt(p[1][0]**2+p[1][1]**2+p[1][2]**2);
P=gRandom->Gaus(pmu,pmu*resoP);
dmy2=10/c/(pmu/sqrt(mmu**2 +pmu**2))*1e9;//ns
dt=gRandom->Gaus(dmy2,dmy2*resoT);
total->Fill(P,dt);}
else if(rndm>Gamma1+Gamma2+Gamma3+Gamma4&&rndm<=Gamma1+Gamma2+Gamma3+Gamma4+Gamma5){
twobodydecay(mK, mpi, ppi, mpi0, ppi0);
P=gRandom->Gaus(ppi,ppi*resoP);
dmy2=10/c/(ppi/sqrt(mpi**2 +ppi**2))*1e9;//ns
dt=gRandom->Gaus(dmy2,dmy2*resoT);
total->Fill(P,dt);
}





32

3.3

3 粒子衝突と崩壊 ∼ 静止 K + 崩壊を中心に ∼静止 K 崩壊荷電粒子のスペクトラム

sample11.cxx の続きelse if(rndm>Gamma1+Gamma2+Gamma3+Gamma4+Gamma5&& rndm<=Gamma1+Gamma2+Gamma3+Gamma4+Gamma5+Gamma6){
DalitzdecayKto3piPlus(mK,mpi,p[0][0],p[0][1],p[0][2],mpi,p[1][0],p[1][1],p[1][2],mpi,p[2][0],p[2][1],p[2][2]);ppi=sqrt(p[0][0]**2+p[0][1]**2+p[0][2]**2);
P=gRandom->Gaus(ppi,ppi*resoP);
dmy2=10/c/(ppi/sqrt(mpi**2 +ppi**2))*1e9;//ns
dt=gRandom->Gaus(dmy2,dmy2*resoT);
total->Fill(P,dt);}
else if(rndm>Gamma1+Gamma2+Gamma3+Gamma4+Gamma5+Gamma6&& rndm<=Gamma1+Gamma2+Gamma3+Gamma4+Gamma5+Gamma6+Gamma7){
DalitzdecayKto3piZero(mK,
mpi0,p[0][0],p[0][1],p[0][2],mpi0,p[1][0],p[1][1],p[1][2],mpi,p[2][0],p[2][1],p[2][2]);ppi=sqrt(p[2][0]**2+p[2][1]**2+p[2][2]**2);
P=gRandom->Gaus(ppi,ppi*resoP);
dmy2=10/c/(ppi/sqrt(mpi**2 +ppi**2))*1e9;//ns
dt=gRandom->Gaus(dmy2,dmy2*resoT);
total->Fill(P,dt);}
else continue;
kaon->Fill();if(i%10000==0)cout<<"."<<flush;}
c1=new TCanvas("c1","",500,300);
total->SetTitle("");
total->SetXTitle("Momenum p (MeV/c)");total->GetXaxis()->CenterTitle();total->SetYTitle("#Delta t (ns)");total->GetYaxis()->CenterTitle();total->Draw("colz");kaon->Write();c1->Print("output.pdf");
return;
}




33

3.3

3 粒子衝突と崩壊 ∼ 静止 K + 崩壊を中心に ∼静止 K 崩壊荷電粒子のスペクトラム

120
110

4

10
100
90

3

10
∆ t (ns)

80
70
2

10

60
50

10

40
30
20
0

50

100

150
Momenum p (MeV/c)

200

250

300

1

図 6 時間分解能 100 ps、運動量分解能 0.9% におけるスペクトラム分布.

図 6 に運動量と時間差のスキャッタ plot を示します。z 軸は Log スケールで表示しています。235 MeV/c あたりに存在する Kµ2 の右にちょこんと存在しているのが Ke2 なのですが、これではわかりませんね。時間分解能を向上∆ t (ns)

すれば見えてくることは明らかで、例えば、20 ps に設定してもう一度見て見ましょう。

5

50

10

45

10

40

10

4

3

2

10

35

10

30

160

図7

180

200
220
Momenum p (MeV/c)

240

260

1

時間分解能 20 ps、運動量分解能 0.9% におけるスペクトラム分布.

図 7 は TOF 時間分解能 20 ps における運動量と時間差の分布です。こうなってくると 245 MeV/c あたりで Ke2が分離しました。さて、時間差閾値をどこに設定するかで Kµ2 によるなだれ込みの量が変わって来ます。ROOT fileにデータを格納しましたので、実際に解析していきましょう。

34

3.3

3 粒子衝突と崩壊 ∼ 静止 K + 崩壊を中心に ∼静止 K 崩壊荷電粒子のスペクトラム

sample12.cxx

void sample12(){gStyle->SetOptStat(0);double P,dt;
double thre=33;
char name[100];
file=new TFile("kaon.root");kaon=(TTree*)file->Get("kaon");kaon->SetBranchAddress("P",&P);kaon->SetBranchAddress("dt",&dt);
TH1F*mom[10];for(int i=0;i<8;i++){sprintf(name,"mom%d",i);mom[i]=new TH1F(name,name,300,0,300);}
for(int i=0;i<kaon->GetEntries();i++){kaon->GetEntry(i);for(int j=0;j<8;j++){if(dt<thre+0.2*j) mom[j]->Fill(P);}if(i%100000==0)cout<<"."<<flush;}
c1=new TCanvas("c1","",500,300);
mom[0]->SetTitle("");
mom[0]->SetXTitle("Momenum p (MeV/c)");mom[0]->GetXaxis()->CenterTitle();mom[0]->GetXaxis()->SetRangeUser(220,260);
mom[0]->GetYaxis()->SetRangeUser(0,50);
mom[0]->Draw();
for(int i=1;i<4;i++){mom[i]->SetLineColor(i+1);mom[i]->Draw("same");
}




35

4

まとめ


sample12.cxx の続きtext=new TLatex(247,20,"K_{e2}");text->SetTextSize(0.05);
text->Draw("same");text=new TLatex(235,45,"K_{#mu2}");text->SetTextSize(0.05);
text->Draw("same");text=new TLatex(225,5,"K_{e3}");text->SetTextSize(0.05);
text->Draw("same");c1->Print("output.pdf");
return;
}




50

Kµ 2

45
40
35
30
25

K e2

20
15
10
5
0
220

K e3
225

図8

230

235

240
245
Momenum p (MeV/c)

250

255

260

TOF 事象選択 ∆t < 33, 33.5, 34, 34.5 ns における運動量分布.

TOF の時間差の事象選択を ∆t < 34.5, 34, 33.5, 33 ns ときつくしていくと、次第に Kµ2 が減っていく。ちょうど ∆t < 33.5 ns の時 (赤)、Kµ2 のなだれ込みを優位に抑制でき、Ke3 と Ke2 のスペクトラムのみが見えるようになりました。ここまですれば Ke2 を観測し、Ke2 の事象数を数えることができます。つまり分岐比の測定ができることを意味します。分岐比 Γ(Kµ2 ) と Γ(Ke2 ) の比を測定することができます。
もしこの値が理論値より統計的優位に離れているなら、暗に新しい物理を示唆します。

4 まとめ
本稿では静止 K による崩壊を中心に模型作成とデータ解析を実施しました。時間分解能 100 ps では、運動量と時間差の分布で Kµ2 のテイルがまだまだ Ke2 に雪崩込んで見えません。20 ps の TOF があれば粒子識別は TOF だけ
で Ke2 を観測することができるが、ちょっと現実的ではありません。実際、Kµ2 のテイルから Ke2 を抜き出すためには、さらにγ線のデータを駆使する必要があります。なぜなら、e+ が標的内もしくは磁場内で制動放射した場合、運動量が下がり Kµ2 の領域になだれ込み、さらに K + → e+ νe γ の３体崩壊が分岐比 9.5 × 10−6 で入り込むため、誤

36

4

まとめ

認する可能性を持つからです。K + → µ+ νµ γ は分岐比 1.33 × 10−5 なのでγ線を放出して TOF で事象選択すれば無視できる水準に抑制できる。e+ の制動放射によるエネルギー損失はどうしようもないので、実際に実験をするときは 105 程度大きい分岐比を持つ Kµ2 を徹底的に排除するために TOF 以外に粒子識別を２重３重にかける必要があ
ります。

あとがき
ROOT 入門 第２弾を書きましたが、第１弾が 2014 年なので 3 年越しなんです。今回は ROOT を勉強するというより、C/C++ を勉強した感が強い気がします。そして何かしらの物理を入れて入門らしくしたいな〜という気持ちもありましたので、どこの実験とは言いませんが、最初のかじりをテーマにしました。執筆していてなんですが、自分にもいい勉強になった気がします。2 体崩壊は簡単なのに 3 体崩壊になると途端に難しくなりました。GEANT4のソースを参考にして Phase Space や Dalitz Decay など組み込んでいます。分岐比は PDG を参考にしています。
それでも物理が間違っているとか、お偉い方々に見られると指摘される箇所があるかもしれません。なので本稿を100% 鵜呑みにしないでください。最後に、本稿で取り扱ったサンプルソースは改造したりして有意義に使ってください（大したことはやっていないのですが）。

37