如何用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 & 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 & projSet, 191const RooAbsData & 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),毕竟是论坛嘛。主要是可以看看大家都在关心什么问题,可以了解一下。
~~继续思考~~ |
-------开放研究之所以可能 |