なぜ加速器を使うのか？〜素粒子と初期宇宙〜①

さて、ここでは加速器や素粒子についておおまかな背景を述べる。

僕は原子核衝突実験の人なので、だいぶそっち方面に偏った説明になるかもしれない。

 

加速器の話

スイス、フランスの国境付近に CERN（欧州原子核研究機構）と呼ばれる研究所があり、せるんとか呼ばれている。CERNにはLHCと呼ばれる山手線一周ほどの巨大な円形加速器がある。

LHCは多くの国と研究機関が力を合わせて１兆円もの大金をかけ作られた大規模なものであり、小さな小さな粒子をほぼ光速度まで加速、衝突させては、世界中の物理学者がキャッキャウフフと喜んでいるのだ。

LHCだけでなく世界中に様々な加速器があり、日本にもJ-PARCやKEKと呼ばれる加速器がある。ここではその種類については詳しくは述べないが、僕が言いたいのは世界中の科学者がこぞって加速器をつくっては高いお金をかけて実験しているということだ。

なぜこんなことをするのだろうか？

一般の人にとっては彼らは気が狂っているとしか思えないだろう。またはお金の無駄遣いにしか思えないだろう。

ここではなぜ彼らがそこまで加速器に熱中するのかを言いたいと思う。

↑CERNの上空写真。加速器は地中に埋まっているので直接は見えないが、とにかくでかい。

 

原子と原子核の話

加速器の話をいったんおこう。

「物質の最小の構成要素はなにか」というのは人類が長い間ずーっと追い求めてきた課題である。誰もが一度は考えたことがあるのではないだろうか？

小学校ですべての物質は原子からできていると習うだろう。

といっても原子も最小の構成要素ではなく、原子はさらに原子核とその周りを回る電子から成る。

原子核はさらに陽子や中性子からできていて、陽子の数で原子の種類が決まる。

例えば水素なら陽子が一個だし、炭素なら１２個である。

中性子の数は関係ないのかって？

まったく関係なくはないのだが、原子の性質はほとんど陽子で決まる。

その理由は中性子は電荷を持たないが、陽子は＋１の電荷を持っていることに起因する。

原子は全体で電荷中性になるのが安定になるため陽子と同じ数だけ電子が外側を回っている。化学反応というのはざっくり言えば電子をやりとりする反応なので、陽子の数=電子の数=原子番号 であり、これが原子の種類を決めることになる。

ここらへんは主に化学の話で、一般の物理を専攻としていない人が学校で習うのはこのへんまでだろう。

↑原子のイメージ、中国電力HPより

素粒子の話

陽子や中性子をまとめて核子と呼ぶが、核子はさらに３つのクォークと呼ばれる素粒子からできており、これが物質の最小の構成要素である。（と現在では考えられている）

この世界に存在するすべての物質は原子からできていることを考えれば、つまるところ物質の最小の構成要素である素粒子の性質がわかればこの世界はどのようにできているのかということを突き止められるかもしれない。

素粒子を研究するということは物質の根源に迫るということなのだ。

素粒子の種類は1つではなくいくつか種類があり、標準模型と呼ばれる理論でまとめられている。

↑標準模型。素粒子は実は１種類じゃなくていろいろある。

レプトンと書いてあるのは電子の仲間で３種類あるニュートリノもレプトンである。e、μ、τの3種のレプトンはどれも-1の電荷を持つが、ニュートリノは電荷も持たなければ質量もほぼ0 であるためほとんどの物質もスカスカすり抜けてしまう。

ほとんどの人にとっては「電子は知ってるけどμ粒子やτ粒子なんて知らねーよ」

って感じだろう。これはμやτは電子に比べて質量が重く不安定なのですぐに崩壊していまうためである。

クォークと書いてあるのが陽子や中性子のような物質を形づくっている素粒子である。右にあるボソンというのは相互作用を媒介する粒子であり、ヒッグスと書いてあるのは物質に質量を与えるヒッグス粒子という割と最近（ようやく）見つかった粒子である。

また、クォークやレプトンにはそれぞれ電荷などの物理量が反転して反粒子というやつがいる。例えば電子の反粒子は陽電子というプラスの電荷を持つもので、電子の陽電子が出会うと対消滅してエネルギーに変わる。我々が普段反粒子を見かけないのはこの世界は粒子の数＞反粒子の数であり、反粒子が生成されてもすぐに対消滅して消えてしまうからだ。SFとかで好まれそうなネタだ。

反粒子が発見の歴史は…...などなどと全部説明しているときりがないのでこれくらいにしておいて、ここでは核子とそれを構成するクォークに注目してみよう。

クォークと強い相互作用の話

陽子や中性子などの核子は3つのクォークから構成されている。

uクォーク２個とdクォーク１個から構成されているのが陽子で、uクォーク１個とdクォーク２個から構成されてるのが中性子である。

上の表を見るとuやd以外にもbとかtとかいろいろあるが、これら前述のμやτ同様質量が大きくすぐに崩壊してしまうので自然界でではなかなかお目にかかれない。

これらの重い粒子は加速器を用いて人工的に創りだすことで実験で初めて発見されたのである。（ただしμ粒子なんかは実は宇宙線として空から降ってきたりしている）

原子核の周りを電子が回っているのは電磁気力によるものだが、クォーク同士や陽子と中性子同士は強い相互作用と呼ばれる力で結びついている。

３つのクォークはグルーオンと呼ばれる素粒子を媒介として結びついている。

媒介として、というのは物理を専門にしていない人にはしっくりこないかもしれない。

例えば電磁気力を例に取ると、プラスの電荷を持った粒子とマイナスの電荷を持った粒子

には引力が働くが、素粒子論だとこれを光子をキャッチボールのようにやり取りすることで引力が働くと考える。光子というのはつまりは光のことである。光は電磁波という波であるが、素粒子の世界ではこれをひとまとまりの粒子のようにとらえ、「光子」と呼んでいるのだ。

同じようにクォークはグルーオンと呼ばれる素粒子を通して結びついているのである。

さて光子と呼ぶとピンとこなくても、電子や光を知らない人はいないだろう。

しかしクォークやグルーオン、強い相互作用と言われても日常生活でそれを感じることはまずない。その理由は強い相互作用が原子核内部のごくごく近距離でしか働かないからである。

原子核にくらべてずっと大きなスケールで生活いている僕らには日々の生活で強い相互作用を感じることはないのだ。

我々が日常生活で感じる「力」は重力を除き、すべて電磁気相互作用によるものだと考えて良い。

また、強い相互作用は力が働く範囲こそ狭いが、電磁気相互作用に比べてめちゃくちゃ強い。名前の由来もそこにある。

放射線や被爆の恐ろしさも実はそのへんに密接に関連している。

我々が普段目にするような化学反応は原子核の周りを回る電子の反応であるが、放射線は原子核の崩壊によって飛び出してくる粒子である。

したがって放射線は通常の化学反応と比べて桁違いのエネルギーを持っているため、人体に浴びてしまうとDNAを傷つけてしまうのである。これが俗に言う被爆だ。

放射線を理解するには原子核の性質を理解することが重要なのである。

ちなみに、中高校生時代僕は陽子はプラスの電荷を持ってて反発し合うはずなのに原子核が結びついているのはおかしくないか？とずっと思っていたが、これも電気的な反発力よりも強い相互作用のほうが文字通り結びつきが「強い」からである。

どんどん脱線してしまうが、この陽子と中性子は中間子と呼ばれる粒子を媒介として結びついており、この中間子の存在を予言してノーベル賞をもらったのがかの有名な湯川秀樹である。

クォークの性質をどうやって調べるか？

原子はとても小さいので顕微鏡で直接みることなんて不可能だった。そのため昔（といっても１００年位前の）人はあれやこれやとがんばって原子の存在を立証したのである。*1ただ最近はからり技術が発達してきて、驚くべきことに原子のスケールでは実際に原子を顕微鏡で見たりつまんで動かしたりできるようになったらしい。IBMという企業が原子の配列を変えて世界最小のアニメーションを作ったのが少し前に話題になったが、これは本当にすごいと思う。

↑原子の配列を変えて作った世界最小のアニメーション

ただ、原子に比べて原子核はさらに小さいのでこんなことはまだできない。

原子も原子核もどっちも小さいのでピンと来ないかもしれないが、原子に比べて原子核はめちゃくちゃ小さい。原子の大きさがだいたいmなのに対して、原子核の大きさはmくらいなので5ケタ違う。よくある例えだと、原子が東京ドームだとすると原子核の大きさはその中に置かれた一円玉くらいの大きさだ。

さらにやっかいなのはクォークの閉じ込めと呼ばれる性質だ。

電磁気の場合はプラスの粒子とマイナスの粒子を遠ざければ遠ざけるほど力は弱まるが、グルーオンの場合クォークの距離を引き伸ばせば引き伸ばすだけ結合が強くなるバネのような性質がある。

クォークを核子から取り出そうとするには無限のエネルギーが必要となってしまい、もはや小さいうんぬんの話ではなく、クォーク１個を取り出して見ることは原理的に不可能なのだ。

これでは八方塞がりであるが、これらを打開する手が実は一つだけあるのである。

つづく。

*1 このへんの原子の発見の歴史はめちゃくちゃ面白い。「だれが原子をみたか」という本がおすすめ。原子の発見の歴史がわかると改めて熱力学や統計力学というものが楽しくなるだろうと思う（僕はなった）。今でこそ熱＝原子の運動というのが常識となっているが、実は熱力学の理論が確立したのは原子の発見よりも前なのである。熱力学が原子の存在を過程しなくてもマクロな物理量だけを使って体系的に成り立っているのはとても興味深い。このへんは田崎晴明さんの「熱力学」「統計力学」なんかがおすすめ。

高エネルギー物理の基礎〜はじめに〜

はじめに

さて、このブログは今までタイトルの通り自分用の備忘録として使ってきた。

PCとかプログラミングの設定などのメモを中心に書いてきたが、せっかくなので自分が専攻している物理学、特に高エネルギー物理学と呼ばれるものについてもまとめておこうと思う。

高エネルギー物理のカテゴリから順番に、教科書のような雰囲気で読めるようにするのが目標である。

なので、順番を変えたり、記事の内容をあとから修正したりということが割と頻繁に起こるかもしれない。

今までの記事は自分用のメモであったが、このカテゴリの記事は多くの人に読んでもらいたいと思っているし、疑問に思うところや間違いなどがあれば積極的に指摘、議論してもらいたいとも思っている。

想定している読者層

高エネルギー物理を専攻している学生だけでなく学部生の物理の学生、素粒子物理学や原子核物理学に興味を持ち始めている物理を専攻していない人にも広く読んでもらいたいと思っている（願望）。

ただ、ブリーバックスやニュートンみたいな基本的に数式が出てこない本と違って、数式を出しながら説明していきたいと考えているので、数式に拒否反応を示す人は読むのが辛くなるかもしれない。数式を読み飛ばしながら読んでもわかるように考慮はするつもりである。

これからの記事をよんでなんとなく「あーこの分野の人たちはこういうことをやっているのかぁ」

ということが伝わればいいと思う。また、このブログは高名な物理学者の人が書いた本と違い、ポンコツの一学生が頑張って書いているものなので、記事の内容を鵜呑みにするのは非常に危険である。

なぜこんな記事を書こうと思ったか

さて、なぜこのようなうぬぼれたシリーズの文章を書こうと思ったのかについてもまとめておこう。僕が大学学部生の時に最も衝撃を受けたサイトが「EMANの物理学」というサイトで、物理を専攻している人なら一度は聞いたことがあるのではないかと思う。

当時大学で使っていた教科書からはどうも堅苦しさを感じていたのだけど、webサイト上でこんなにもわかりやすく、面白く物理を語ることができるのかと当時の僕は夢中になり、サイト内の文章はすべて読破した。今の自分が物理を好きだと言えるのはかなりこのサイトの影響があると思う。それと同時に自分が勉強した物理について文章を書いてもこのサイトの二番煎じや下位互換にしかならないと思い、自分から積極的にこのような文章を書こうとは思いもしなかった。

きっかけは修士論文

さて、満を持して修士課程に進んだ僕であったが、そこでは高エネルギー物理の分野で初めて触れるような式や用語がたくさんあった。そして、それらに関連する文章を読んで「わかり辛いなぁ」と思うことが多々あった。

先に言い訳をしておくと、別にその分野の本や論文を書いている人の説明する力がないわけでなく、専門家の脳内ではちゃんと理路整然と理解できているのだと思うが、頭の弱い僕にはそれがどうにもわかりづらい文章に感じてしまったのである。

自分のような馬鹿や、普段論文を頻繁に読んでいない層にもわかるような文章を書けならいいな、と思う。そしてそれをきっかけとして少しでも物理に興味を持ってもらえれば、とも思う。

修士論文を書くときに指導教官の一人から「修士論文は未来の後輩に向けた教科書を書くような感覚で書くといいよ」というようなことを言われて張り切ったのだけど、仮にも「論文」と名前のつく性質上あまりいい加減は文章はかけず、どうにももやもやした気持ちが残った。

このブログはそんな自分の憂さ晴らしである。

偉そうなことをいっぱい書いたが、まだまだ模索している段階なので暖かく見守ってもらえればと思う。

丸め誤差の恐怖〜浮動小数点数の扱いに注意する〜

丸め誤差は予想以上に厄介であることを思い知らされた。

こいつのせいで予想外のバグとの格闘に時間を割かれてしまったのだ。

自分への戒めとしてこの文章を残す。

 

丸め誤差とは？

簡単な例として10/3を考えよう。10を3で割り切ると

3.33333....

と永遠に続き割り切れない。これらの情報すべてを正確に残すには無限の容量が必要なので適当なところで打ち切る必要がある。有効数字10ケタまでを情報として残そう！みたいな。

こうやって有効数字の桁数で切り捨てた時に発生する誤差が丸め誤差である。

ただ、これは10進数の場合であり、コンピュータは2進数なのでちょっと異なる。

例えば0.１は10進数だと誤差なしで書けるけど、２進数で表現すると0.000110011...みたいな循環少数になってしまう。

僕達はふだんついつい10進数脳で物事を考えてしまう上に、丸め誤差もあまり気にせずプログラムを書きがちだけど、これは要注意である。

例

自分がハマった例を書く。

x、yと2つの量(3つの要素をもつベクトル)があったとして

ｘ=(3,3,4)

ｙ=(1,2,7)

とする。ここでx-yの期待値を考えてみよう。ｘもyも3つの要素をすべて足すと10になるのでx-yの期待値は0になるはずである。実際に手計算してみると

{(3-1)+(3-2)+(4-7)}/3 = 0

まぁ当たり前の結果だ。

次にこれをc言語で書いてみる。

double a = 3/3. + 3/3. + 4/3.;

double b = 1/3. + 2/3. + 7/3.;

cout<<a-b<<endl;

んで出力は0になるかというと思えば

-4.440892e-16

とか出てくるのである。0じゃない。

３で割った時に丸め誤差がab間で発生してしまい、差をとっても0にならなくなってしまうのだ。

「10の-16乗でしょ？ほとんど０なんだから問題ないじゃん」

とか思うかもしれないがそうでもない。

例えばa-bの逆数をとって

c=1/(a-b)

を計算したいとする。

(a-b)が0だとこの量は発散してしまうので、

if((a-b)!=0) c = 1./(a-b);

とか書けばa-bが０じゃないときだけ計算してくれることになるが、先ほどの値を使うと0じゃないので

普通にこのif文をtrueとみなしてしまい、結果とんでもない値がcに代入されてしまう。

解決策

足してから割ればうまくいった。

double a = (3 + 3 + 4)/3.;

double b = (1 + 2 + 7)/3.;

cout<<a-b<<endl;

一度足してしまえば割り算は一回なので、誤差が発生するリスクも少ない。

今回の場合はどちらも10/3.を計算することになるのでぴったり一致する。

誤差の桁の意味

せっかくなのでもう少し踏み込んでみよう。

最初の例だと

a-b=-4.440892e-16

なので、15桁は一致していて、16桁目がずれていることになる。

今回は浮動小数点数のdoubleを使っているが、doubleの場合

（符号1ビット）*(仮数部52ビット)*(指数部11ビット)

の62bit = 8byteからなる。

今回の例だと指数部は問題ないので、大事なのは仮数部の部分で、ここが有効数字の桁数を決めてると考えて良い。

2^52＝4.5035996e+15

なので、仮数部は１０進数で15桁くらいの精度で計算できることになる。

そう考えると16桁のオーダーで誤差が発生するのもうなずける気がする。

んーと、合ってるかな？？

＊余談

丸め誤差と似たような量として打ち切り誤差ってのもある。

これサイトによって書いてあることが違って、あるサイ卜だと

「無限級数などの無限に続く計算を途中で打ち切った時に生じる誤差」

っとか書いてある一方で

「循環少数などの計算を途中で打ち切った時に生じる誤差」

って書いてあるページもある。

前者だと足し算の項とかを足すのを途中でやめてるわけだから今回の例とはちょっと異なるし、後者だと今回の無限小数の例も打ち切り誤差ってことになる。

どっちなんだ。手元の本だと前者に近くて、丸め誤差と打ち切り誤差は区別されている。

まぁ無限に続いているものを途中で打ち切ることで生じる誤差、という意味では丸め誤差と似たような量なんだろう。

他にも情報落ち、ケタ落ち、オーバーフローとかいろいろあるけど、今回はこのへんで。

まとめ

・割ってから足すな、足してから割れ

おわりー

g++でc++11の機能を使ってコンパイルする時のメモ

コンパイルの方法

g++ -std=c++11 -O2 -Wall ‘root-config --cflags —libs’ -o test test.C

とすればtestという名前のオブジェクトファイルができる。

./test

とかすれば実行される。

main関数に引数を与える場合は

./test (引数1)  (引数2)　とかやれば通る。

一度書いてしまえば毎回コピペなのだけど、一応オプションの意味も書いておく。

基本形は

g++ -o test test.C

でtest.Cをコンパイルして、testという名前という名前の実行ファイルを作ってくれる。

-O2 : 最適化を行なってくれる。最適化を行うと、コンパイルの時間はかかるが実行時間とかメモリ使用量が軽減されるらしい。

-Wall : 警告を表示するオプション

‘root-config —flags —libs : rootを使う時に必要。（rootはCERNが開発しているデータ解析フレームワークです）

-std=c++11 : c++11の機能を有効にする

コンパイラのバージョンの確認

unorderd mapなど、c++11から実装された新しいテンプレートクラスを使うにはコンパイラもそれに対応していなければならない。

環境によってコンパイルが通らないことがある。

$ gcc -dumversion

でバージョンの確認。

4.8.2

とか出てくる。もうちょい詳しく見たいときは

$gcc —version

とすれば詳細がいろいろ出てくる。

gccの部分をg++にすればg++のバージョンも出てくる。

どちらもほぼ同じだった。

このg++は2014年と書いてある。この場合はc++11の機能は使えた。

一方で4.4.7 (201203013)のほうだとこけた。

cc1plus: error: unrecognized command line option “-std=c++11”

とでる。「c++11なんて知らんわ」というエラー。

https://gcc.gnu.org/projects/cxx-status.html#cxx11

このページを見た限り、バージョン4.7だか4.8以降じゃないと使えないみたいだ。

以上！

c++ でテキストファイルを読み込んでループを回す時のメモ

データ解析においてテキストファイルに書かれたデータを読み込む機会は結構多い。

俺も今までさんざん使ってきたんだけど、結構細かい書き方忘れちゃうし、曖昧な箇所も多いのでこの際ちゃんと覚える。

例えば次のような内容のテキストファイル（test.dat）があったとする。

1  2

3  4

これを読み込んでみよう。

書き方その1

 void test(){
    int a,b;
    ifstream ifs;
    ifs.open("test.dat");
    int loop=0;
    ifs>>a>>b;
    while(!ifs.eof()){

    cout<<"loop="<<loop<<endl;
    cout<<"a,b="<<a<<","<<b<<endl;
    ifs>>a>>b;
    loop++;
  }

}

出力結果

loop=0

a,b=1,2

loop=1

a,b=3,4

メモ

ヘッダはfstream.hが必要。using name spaceとかは書いてないけどコンパイルするならほんとは必要。ファイルの読み込みはfstreamを通して行う。

ファイルを開くとき、ここでは一度fstreamクラスのオブジェクト（ifs）をつくってから、openという関数を使っているが、ifstream("test.dat")と一気にコンストラクタで開くことも可能。まぁどっちでもいい。

loopはloopの回数数えてるだけので別になくてもいい。ただ、for分と違って無限loopが怖いのでwhile文の場合はloop回数についつい神経質になってしまう。

eof()は読みこむ際に最後の行であるかを判定してtrue（０）かfalse（１）を返す。

うまく読み込めればfalse、最後までいって読み込めなくなったらtrueを返す。

  while文はtrueの時だけ回るので!がついているというわけ。

コードをよくよく見ると、１行目をloopの外で読み込んでからloopに突入し、表示→読み込みの順番でloopを繰り返しているが、これにはちゃんと意味がある。下の悪い例を見てほしい。

悪い書き方

 void test(){
    int a,b;
    ifstream ifs;
    ifs.open("test.dat");
    int loop=0;
    while(!ifs.eof()){

   ifs>>a>>b;

    cout<<"loop="<<loop<<endl;
    cout<<"a,b="<<a<<","<<b<<endl;
    loop++;
  }

}

出力結果

loop=0

a,b=1,2

loop=1

a,b=3,4

loop=2

a,b=3,4

違いはloopの中で読み込み→表示を繰り返してること。見てわかるように最後に一回余分にloopが回ってる。

この不具合は!eof()が最後の行を読み込んだときはtrueでその次の存在しない行を読み込んだ時に初めてfalseを出力するために発生する。

つまり、最後の存在しない行を読み込んで!eof()がfalseになった直後に表示を余分にしてしまうのだ。このへんはなんとなくコードを書いているとやってしまいがちなので注意する。（というかやってしまった）

書き方その２

 void test(){
    int a,b;
    ifstream ifs;
    ifs.open("test.dat");
    int loop=0;
    while(ifs>>a>>b){

    cout<<"loop="<<loop<<endl;
    cout<<"a,b="<<a<<","<<b<<endl;
    loop++;
  }

}

違いはwhile文の中身で、ifs>>a>>bの読み込みがうまくいってる限りloopが回り、読み込めなくなった時点で止まる。さっきよりすっきりしてるし、例の一回余分に読み込んでしまう不具合も発生しない。難点を挙げるとすれば読み込みのifs>>...の部分が長くなったときにwhile文のかっこの中身がごちゃごちゃしてしまって見づらいくらい？いやほんとに難点って言っていいのかわからないけど。

ちなみに言うまでもないけど、for文はloopの回数がわかってる時には使えるが、

datファイルの中身が何行あるかわからないときは今回みたいにwhile文を使うしかない。

個人的にはfor文のほうが好きだけど致し方ない。

まとめ

・テキストファイルの読み込みはfstreamを使いwhile文でloopを回す

・fstream::eof()を使うときは最後の行の振る舞いに注意

おわり！

Linuxで自分の使ってるシェルの種類を確認するときのメモ

自分の使ってるシェルを確認

$ echo $SHELL

で現在使ってるシェルの種類がわかる。自分の環境だと。

/bin/tcsh

とでた。

cshを拡張、バグ修正したものがtcshらしい。

csh、tcshはcsh系と呼ばれ、

sh、bash、kshはsh系と呼ばれる。文法とかも微妙に違う。

zshはどっちもいけるらしい。

$ cat /etc/shells　で使えるシェルを確認。自分の環境だと

/bin/sh

/bin/bash

/sbin/nologin

/bin/ksh

/bin/zsh

/bin/tcsh

/bin/csh

/usr/bin/ksh

/usr/local/bin/tcsh

と出てくる。大抵のシェルは使えるみたいだ。cshの文法で適当にスクリプト書いてみる

test.csh

#! /bin/csh

set a=0

echo $a

csh ./test.csh

で0が出力される。bashでやりたければ

a=0

echo $a

で

sh ./test.csh　とか

bash ./test.csh

とかやればbashでもまわる。拡張子が.cshでも中身がbashの文法であれば普通にbashとして

まわった。どゆこと？と思ったが、どうやら拡張子は見やすくするために付けるだけであまり意味を持たないらしい。

まぁでも拡張子が.cshなのに中身がbashなのは混乱を招くので、test.shとかしといたほうがいいな。

ちなみに最初の#! /bin/cshもcshで回す時に必要だと思ってたけど、これも別になくてもいいらしい。

次に実行権限を与える。

$ chmod +x test.csh

これで

./test.csh

だけで実行できるようになる。#! bin/cｓｈでcshのpathを指定しているのでcshで実行できるみたいだ。

今まで理解がなあなあで使ってた部分がちょっとすっきりした。

rootのTProfileで任意の幅でrebinしたりProfileどうしを足したり割ったりする方法

いくつかのTProfile（or histogram）同士でrebinしたりProfile同士を足したり割ったりしたい。

今まではGetBinContent()でビンの中身持って来て足したりとか脳筋っぽいことをしてたんだけど、

histogramと違ってProfileの場合はエラーの計算とかがめんどくさい。

ヒストはエントリー数足してルートとればエラーになるけど、profileは誤差伝搬とか面倒だからね。

そんなわけで便利な関数がないか調べてみた。

１、TProfile同士を足す。

http://www-he.scphys.kyoto-u.ac.jp/member/n.kamo/wiki/doku.php?id=study:software:root:main

このページに書いてあった。

hist->Add(hist1,x);//各binの値=現在の値+x*hist1の値

重みとか気にしないならx=1がデフォルトなので引数一個だけ入れれば大丈夫。

例えば BeforeProfile[N] というTProfileのオブジェクトを10個足したい場合

TProfile *profile;

for(int i=0;i<N;i++){

  if(i==0)profile = BeforeProfile[i];

  else profile->Add(BeforeProfile[i]);

}

みたいな感じで書けばよい。最初のTProfileには何も入ってないのでfor文の一回目は普通に代入して

2回目からAddしてるのがミソ。最初にこれを忘れて見事にbreakした。初期化って大事だね。

２、Rebinする

次にRebinする。今回は等間隔じゃなくて自分で任意の幅を決めてrebinしたい。

これもさっきのページに書いてあった。

・等分割の場合

TH1D h_rebin* = h1->Rebin(n);         //(1)

TH1D h_rebin* = h1->Rebin(n,"h_rebin"); //第二引数は新しく作るhistoの名前

・好きな幅で分割する

double xbins[n+1] = {};

TH1D h_rebin* = h1->Rebin(n,"h_rebin",xbins);

第3引数に配列をぶっこむことですきな場所でrebinできるらしい。

(1)はn個で1つのbinにするという意味。全体でn個のbinに分けるという意味ではないので注意。

逆に(2)のnは何個に分けるかという意味なのでややこしい。配列の中身がn+1になってるのは最後の値の上限値があるから。

つまり0-100を2分割したかったら{0,50,100}と3つ値が必要ということ。

(2)の場合、TProfileをRebinしても出力はTH1の形式になるので注意。

３、TProfile同士で割る

これもさっきのAddの代わりにDivide()という関数を使ったら値は出たのだが、エラーの値が超でかい。。

よくみるとWarningが出ていて

"WARNING!!: The algorithm in TProfile::Divide computing the errors is not accurate

  Instead of Divide(TProfile *h1, TProfile *h2, do:

  TH1D *p1 = h1->ProjectionX();

  TH1D *p2 = h2->ProjectionX();

   p1->Divide(p2);

とか出てる。Profile::Divideだと正しくエラー計算できないってこと？どういうこっちゃ。

もともとのソースコード見るとコメントに

//THE ALGORITHM COMPUTING THE ERRORS IS WRONG. HELP REQUIRED

help requiredって…

どうやらDivide()は使い物にならなそう。GetBinContentとかでビンの値とエラーを持って来て

誤差伝搬とかで計算するしかなさそう。

結論、TProfileでもAddしたりRebinしたりする関数はあるが、Divideは使いものにならない。掛け算とかする関数もあるみたいだが、使ってないので使えるかわかんない。