My Avatar

Xiaofei

信仰共产主义,后现代主义者,结构主义者,奇妙发现世界~~
`

如何用ROOT做联合拟合?

2019年05月30日 星期一, 发表于 北京 | 浏览量: -

如果你对本文有任何的建议或疑问, 谢谢! :)

“你本来不该知道这些的,可是现在我们没有更好的工具,所以得这样弄一下……”

如何做联合拟合?

首先这里就是应用Roofit的例子(Tutorials)来学习做联合拟合(SimutaneousFit),这是一个比较常用的东西,当利用多个衰变道来重建共同的共振态时,可能要share一些变量,这样在拟合中可以得到描述该共振态共同的变量值,进而完成拟合。

这里就试着来看看例子

  1首先这里是个声明,在总的ROOT体系里使用了[LGPL2.1协议](https://root.cern.ch/license),这位WouterVerkerke贡献了该例子
  2//////////////////////////////////////////////////////////////////////////
  3//
  4// 'ORGANIZATION AND SIMULTANEOUS FITS' RooFit tutorial macro #501
  5//
  6// Using simultaneous p.d.f.s to describe simultaneous fits to multiple
  7// datasets
  8//
  9// 07/2008 - Wouter Verkerke
 10//
 11/////////////////////////////////////////////////////////////////////////
 12
 13// 头文件,名字起的很有代表性,大家用的程序要有编程规范,不然容易乱套
 14// 其实主要是需要编译时最需要,不过编译好久没有用了,现在我不太会了忘记了,一般找个MakeFile比较好
 15// CMake是用来写MakeFile的编译工具,很有意思,MakeFile用来指挥编译按步骤进行,CMake用来写指挥MakeFile如何写
 16
 17#ifndef __CINT__
 18#include "RooGlobalFunc.h"
 19#endif
 20#include "RooRealVar.h"
 21#include "RooDataSet.h"
 22#include "RooGaussian.h"
 23#include "RooConstVar.h"
 24#include "RooChebychev.h"
 25#include "RooAddPdf.h"
 26#include "RooSimultaneous.h"
 27#include "RooCategory.h"
 28#include "TCanvas.h"
 29#include "TAxis.h"
 30#include "RooPlot.h"
 31
 32// 命名空间,像一种局部用的目录,防止重名,当然没遇到重的所以暂时还不懂
 33using namespace RooFit ;
 34
 35void rf501_simultaneouspdf()
 36{
 37  // C r e a t e   m o d e l   f o r   p h y s i c s   s a m p l e
 38  // -------------------------------------------------------------
 39  // Create observables
 40  // 这里首先定义了变量x用来构建实验分布的模型
 41  RooRealVar x("x","x",-8,8) ;
 42
 43  // Construct signal pdf
 44  // 定义了信号的概率密度函数,这里有中心值、sigma自变量x一起构造高斯函数
 45  // 经常有听到拿“形状”,指的其实就是概率密度函数的形状
 46  RooRealVar mean("mean","mean",0,-8,8) ;
 47  RooRealVar sigma("sigma","sigma",0.3,0.1,10) ;
 48  RooGaussian gx("gx","gx",x,mean,sigma) ;
 49
 50  // Construct background pdf
 51  // 类似于前面构造二阶切比雪夫作为本底概率密度函数,a0和a1作为其中系数
 52  RooRealVar a0("a0","a0",-0.1,-1,1) ;
 53  RooRealVar a1("a1","a1",0.004,-1,1) ;
 54  RooChebychev px("px","px",x,RooArgSet(a0,a1)) ;
 55
 56  // Construct composite pdf
 57  // 再次定义变量f,作为最后构造的物理分布模型的自变量,所得值为上两函数的和,用了RooArgList类
 58  RooRealVar f("f","f",0.2,0.,1.);
 59  RooAddPdf model("model","model",RooArgList(gx,px),f);
 60
 61  // C r e a t e   m o d e l   f o r   c o n t r o l   s a m p l e
 62  // --------------------------------------------------------------
 63  // Construct signal pdf.
 64  // NOTE that sigma is shared with the signal sample model
 65  RooRealVar mean_ctl("mean_ctl","mean_ctl",-3,-8,8) ;
 66  RooGaussian gx_ctl("gx_ctl","gx_ctl",x,mean_ctl,sigma) ;
 67
 68  // Construct the background pdf
 69  RooRealVar a0_ctl("a0_ctl","a0_ctl",-0.1,-1,1) ;
 70  RooRealVar a1_ctl("a1_ctl","a1_ctl",0.5,-0.1,1) ;
 71  RooChebychev px_ctl("px_ctl","px_ctl",x,RooArgSet(a0_ctl,a1_ctl)) ;
 72
 73  // Construct the composite model
 74  // 这里可以看到,不同的“形状”代表了不同的实验分布,分别反映了信号概率密度函数形状,探测器分辨等
 75  // 然后将他们“组合”也就是卷积在一起,来表达该物理过程我们想要了解的实验分布
 76  RooRealVar f_ctl("f_ctl","f_ctl",0.5,0.,1.) ;
 77  RooAddPdf model_ctl("model_ctl","model_ctl",RooArgList(gx_ctl,px_ctl),f_ctl) ;
 78
 79  // G e n e r a t e   e v e n t s   f o r   b o t h   s a m p l e s
 80  // ---------------------------------------------------------------
 81  // Generate 1000 events in x and y from model
 82  RooDataSet *data = model.generate(RooArgSet(x),100) ;
 83  RooDataSet *data_ctl = model_ctl.generate(RooArgSet(x),2000) ;
 84
 85  // C r e a t e   i n d e x   c a t e g o r y   a n d   j o i n   s a m p l e s
 86  // ---------------------------------------------------------------------------
 87  // Define category to distinguish physics and control samples events
 88  // 这里定义了一个目录其中命名均为字符串,所以名字里面空格会识别为字符,所以要谨慎,加不加空格是两个名字
 89  RooCategory sample("sample","sample") ;
 90  sample.defineType("physics") ;
 91  sample.defineType("control") ;
 92
 93  // Construct combined dataset in (x,sample)
 94  // 将数据导入到RooDataSet里面,其中有目录Index(sample),还有变量x,并将目录中的名字与数据集关联起来
 95  RooDataSet combData("combData","combined data",x,Index(sample),Import("physics",*data),Import("control",*data_ctl)) ;
 96
 97  // C o n s t r u c t   a   s i m u l t a n e o u s   p d f   i n   ( x , s a m p l e )
 98  // -----------------------------------------------------------------------------------
 99  // Construct a simultaneous pdf using category sample as index
100  // 这里我们将之前定义的目录送入RooSimultaneous类中去,然后
101  RooSimultaneous simPdf("simPdf","simultaneous pdf",sample) ;
102
103  // 这里将之前定义的分布模型放到RooSimultaneous类的sample里面去,
104  // Associate model with the physics state and model_ctl with the control state
105  simPdf.addPdf(model,"physics") ;
106  simPdf.addPdf(model_ctl,"control") ;
107
108  // 这里将数据进行拟合,拟合的目标物理分布的模型就是前面的
109  // P e r f o r m   a   s i m u l t a n e o u s   f i t
110  // ---------------------------------------------------
111  // Perform simultaneous fit of model to data and model_ctl to data_ctl
112  simPdf.fitTo(combData) ;
113  // 到此为止,就完成了拟合的计算过程,其实已经结束了,只是不直观要画图
114  // 可以看出,RooSimultaneous就是一个类中包含了多个模型,这些多个模型的公共参数用来描绘共同的一个物理分布。
115  // 这种多个模型,可以来自于不同衰变道,也可以来自于不同的样本如例子中随机产生的物理分布和控制样本。
116
117  // P l o t   m o d e l   s l i c e s   o n   d a t a    s l i c e s
118  // ----------------------------------------------------------------
119  // Make a frame for the physics sample
120  // 画图的RooPlot,用来载入画图的内容,在这里定义了分bin
121  RooPlot* frame1 = x.frame(Bins(30),Title("Physics sample")) ;
122
123  // Plot all data tagged as physics sample
124  // 将数据画在上面定义的frame1中,并且通过Cut条件只看physics目录的内容
125  combData.plotOn(frame1,Cut("sample==sample::physics")) ;
126
127  // Plot "physics" slice of simultaneous pdf.
128  // NBL You _must_ project the sample index category with data using ProjWData
129  // as a RooSimultaneous makes no prediction on the shape in the index category
130  // and can thus not be integrated
131  simPdf.plotOn(frame1,Slice(sample,"physics"),ProjWData(sample,combData)) ;
132  // 这里可以看到是又画了一次,这次是单独画的本底用的kDashed的线,可以按后面给的方法搜索类型继续看
133  simPdf.plotOn(frame1,Slice(sample,"physics"),Components("px"),ProjWData(sample,combData),LineStyle(kDashed)) ;
134
135  // 定义第二个图的画图的frame
136  // The same plot for the control sample slice
137  RooPlot* frame2 = x.frame(Bins(30),Title("Control sample")) ;
138  combData.plotOn(frame2,Cut("sample==sample::control")) ;
139  simPdf.plotOn(frame2,Slice(sample,"control"),ProjWData(sample,combData)) ;
140  simPdf.plotOn(frame2,Slice(sample,"control"),Components("px_ctl"),ProjWData(sample,combData),LineStyle(kDashed)) ;
141
142  // 定义好画布,将前面的frame调用Draw()成员函数,就画图了
143  TCanvas* c = new TCanvas("rf501_simultaneouspdf","rf403_simultaneouspdf",800,400);
144  c->Divide(2) ;
145  c->cd(1) ; gPad->SetLeftMargin(0.15) ; frame1->GetYaxis()->SetTitleOffset(1.4) ; frame1->Draw() ;
146  c->cd(2) ; gPad->SetLeftMargin(0.15) ; frame2->GetYaxis()->SetTitleOffset(1.4) ; frame2->Draw() ;
147}
148`</pre>
149
150实际上,我们看完上面的内容就可以了解到**其实我们什么都没有做!**我们做的仅仅是:
151
152<pre>`头文件
153函数{
154  定义预期的多个模型model(RooAddPdf)
155  数据data(DataSet),按照样本目录加载
156  调用联合拟合类(RooSimultaneous).fitTo(data)
157  画图,信号本底
158}
159`</pre>
160
161其他的当然也是类似,是计算机帮我们完成了这些操作,当然如涉及到其背后是如何工作的则要继续了解simultaneous背后的数学是如何工作的,比如具体的误差它是如何计算的?整体上参数是如何传递,计算是如何迭代的?编程?不存在的,可以放后去了解。这也学就是我们称之为脚本的原因吧,之所以使用脚本文件,就是相应的程序包提供了简易的借口,可以对应用的用户功能进行自定义。那么有没有办法不写这么多字呢好头大啊?有的,有的人用Python来自动写脚本,有的人直接用PyROOT。
162
163**使用Vim**
164
165如何使用vim?
166
167这个实际上吧vim的手册简单看一下,然后了解一下计算机发展过程,“文本”这个概念是如何发展的,就可以更好地使用vim,我自己在两年用了杨振伟老师推荐的一个简单的vimrc文件[易水博客](http://easwy.com/),还有重新编译了vim用来安装了[vim-snippets](https://github.com/cuixiaofei/vim-snippets)后,再没有怎么搞过了,我想snippets应该已经功能很强大了。
168想起来第一次用vim其实就和我小时候玩游戏一样,上来就键盘挨个按一遍,软件也是,把菜单打开所有按钮按一遍,这就是我从来看任何计算机相关东西的所谓“窍门”。
169
170**继续试着学ROOT有什么方便的方法?**
171
172如果有机会找人教,可以问一下别人,或者有人给脚本,就可以学的很快,因为大家公用的东西都是比较成熟的。在简单用了之后,还要继续学懂一些,或者看到别的ROOT脚本文件有时候发现看不懂,甚至受到了惊吓,下面的一点内容我们试着想一想办法。
173
174首先想到的就是,看一看例子tutorial,ROOT官方带有许多的例子,只要稍加更改就可以学到很多东西,上面篇幅中做联合拟合的内容就是依据这样的方式来做的,下面还有。
175
176一个是搜索使用的类的定义和继承关系,在网站这里[ClassIndex](https://root.cern.ch/doc/master/annotated.html),比如可以试着搜索ProjWData这个成员,也可以在[ClassList](https://root.cern.ch/doc/master/annotated.html)用搜索框或者是按类找名字ctrl+f,这样搜索可以找到ROOT编程的详细内容。
177可以看到,其定义有两个如下:
178
179第一个定义可以看出,其可有两个变量,一是数据集RooDataSet(RooAbsData三个成员中的一个),从下面的继承图可以清楚地了解到程序的调用关系,第二个可能用于分bin的问题,![](/assets/img/RooAbsData.png)
180
181<pre>`RooFit::ProjWData  ( const RooAbsData &amp; projData,
182Bool_t  binData = kFALSE
183)
184`</pre>
185
186从这里上下两个定义可以看出,ROOT软件包整体内部命名是十分规范的,我们在使用的时候也要尽量把变量函数名写的清楚易懂有一定的规范。如BES合作组的BOSS分析环境,可以查看BOSS定义的基本的类和继承图[BOSS Class Reference Documentation](http://bes3.to.infn.it/Boss/7.0.2/html/classes.html),这里可以清楚看到很多分析中使用的内容。当然,如果有较为统一的规范是不容易的,因为在传统研究中各种不同群体都有可能产生各自内部的“独特语言系统”。
187
188第二个定义可以看出,其可有三个变量,一是从数据集里取数据的方法(RooCategory是其中一个),第二个见上,第三个可能用于分bin的问题
189
190<pre>`RooFit::ProjWData  ( const RooArgSet &amp; projSet,
191const RooAbsData &amp; projData,
192Bool_t binData = kFALSE
193)
194`</pre>
195
196可以看出我们这里使用了两次,一次第一个定义,一次第二个定义,并且没有用bin变量。
197这样我们通过查看这些的定义关系,有助于梳理清楚他们的调用关系,用起来就不会发怵。
198所谓计算机程序就是一大堆的相互调用,从ifelse编译一直调用到逻辑门基本不会出错,如果出现问题,就查看定义看哪里调用跑丢了,如果能够解决好把调用的类型放对了,一般程序都能正常运行,就像喂小猫一样。
199如果想要训练一下编程时候的脑子,其实有更直观一点的编程比如说用G语言LabVIEW,你可以看数据流,并且简单的编程结构也比较直观,实际应用也比较广泛,比如BESIII值班室的用户界面系统,可以看到很多都是使用LabVIEW编的。另外甚至可以通过一些软件也能训练,如[Zachronic](http://www.zachtronics.com)公司旗下的软件都很硬核,还有辅助学习各种编程语言的网站[酷代码](www.codingame.com)非常有趣,很全各种计算机语言都能找到。
200
201可以在网上[ROOT Forum](https://root-forum.cern.ch)看别人遇到了什么问题,如果身边没人也可以试着问,还可以看看同行都在关心什么问题?
202这里提供一个问问题模板[参考](https://root-forum.cern.ch/t/how-to-automatically-generate-classes-code/4667):
203
204<pre>`Hi,
205I want to do sth.
206I have learned "sth. or some class" which helps me in doing sth. This looks like it could help me but I don't know HOW.
207I need to to do similar processing but for a different "sth.? dataset?".
208Could some body guide me to the Documentation or Tutorials for doing this?
209
210Thanks
211"yourname"
212
213Sample: "your sample idea"
214Code:"your sample code"

可以找一个改改用,也可以直接随便问(English),毕竟是论坛嘛。主要是可以看看大家都在关心什么问题,可以了解一下。

~~继续思考~~

-------开放研究之所以可能