My Avatar

Xiaofei

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

如何用ROOT做联合拟合?

2019年05月30日 星期一, 发表于 北京

如果你对本小站有任何的建议或者疑问, 可以在 提交Issues, 谢谢! :)

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

如何做联合拟合?

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

这里就试着来看看例子

首先这里是个声明,在总的ROOT体系里使用了[LGPL2.1协议](https://root.cern.ch/license),这位WouterVerkerke贡献了该例子
//////////////////////////////////////////////////////////////////////////
//
// 'ORGANIZATION AND SIMULTANEOUS FITS' RooFit tutorial macro #501
//
// Using simultaneous p.d.f.s to describe simultaneous fits to multiple
// datasets
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////

// 头文件,名字起的很有代表性,大家用的程序要有编程规范,不然容易乱套
// 其实主要是需要编译时最需要,不过编译好久没有用了,现在我不太会了忘记了,一般找个MakeFile比较好
// CMake是用来写MakeFile的编译工具,很有意思,MakeFile用来指挥编译按步骤进行,CMake用来写指挥MakeFile如何写

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooConstVar.h"
#include "RooChebychev.h"
#include "RooAddPdf.h"
#include "RooSimultaneous.h"
#include "RooCategory.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "RooPlot.h"

// 命名空间,像一种局部用的目录,防止重名,当然没遇到重的所以暂时还不懂
using namespace RooFit ;

void rf501_simultaneouspdf()
{
  // 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
  // -------------------------------------------------------------
  // Create observables
  // 这里首先定义了变量x用来构建实验分布的模型
  RooRealVar x("x","x",-8,8) ;

  // Construct signal pdf
  // 定义了信号的概率密度函数,这里有中心值、sigma自变量x一起构造高斯函数
  // 经常有听到拿“形状”,指的其实就是概率密度函数的形状
  RooRealVar mean("mean","mean",0,-8,8) ;
  RooRealVar sigma("sigma","sigma",0.3,0.1,10) ;
  RooGaussian gx("gx","gx",x,mean,sigma) ;

  // Construct background pdf
  // 类似于前面构造二阶切比雪夫作为本底概率密度函数,a0和a1作为其中系数
  RooRealVar a0("a0","a0",-0.1,-1,1) ;
  RooRealVar a1("a1","a1",0.004,-1,1) ;
  RooChebychev px("px","px",x,RooArgSet(a0,a1)) ;

  // Construct composite pdf
  // 再次定义变量f,作为最后构造的物理分布模型的自变量,所得值为上两函数的和,用了RooArgList类
  RooRealVar f("f","f",0.2,0.,1.);
  RooAddPdf model("model","model",RooArgList(gx,px),f);

  // 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
  // --------------------------------------------------------------
  // Construct signal pdf.
  // NOTE that sigma is shared with the signal sample model
  RooRealVar mean_ctl("mean_ctl","mean_ctl",-3,-8,8) ;
  RooGaussian gx_ctl("gx_ctl","gx_ctl",x,mean_ctl,sigma) ;

  // Construct the background pdf
  RooRealVar a0_ctl("a0_ctl","a0_ctl",-0.1,-1,1) ;
  RooRealVar a1_ctl("a1_ctl","a1_ctl",0.5,-0.1,1) ;
  RooChebychev px_ctl("px_ctl","px_ctl",x,RooArgSet(a0_ctl,a1_ctl)) ;

  // Construct the composite model
  // 这里可以看到,不同的“形状”代表了不同的实验分布,分别反映了信号概率密度函数形状,探测器分辨等
  // 然后将他们“组合”也就是卷积在一起,来表达该物理过程我们想要了解的实验分布
  RooRealVar f_ctl("f_ctl","f_ctl",0.5,0.,1.) ;
  RooAddPdf model_ctl("model_ctl","model_ctl",RooArgList(gx_ctl,px_ctl),f_ctl) ;

  // 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
  // ---------------------------------------------------------------
  // Generate 1000 events in x and y from model
  RooDataSet *data = model.generate(RooArgSet(x),100) ;
  RooDataSet *data_ctl = model_ctl.generate(RooArgSet(x),2000) ;

  // 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
  // ---------------------------------------------------------------------------
  // Define category to distinguish physics and control samples events
  // 这里定义了一个目录其中命名均为字符串,所以名字里面空格会识别为字符,所以要谨慎,加不加空格是两个名字
  RooCategory sample("sample","sample") ;
  sample.defineType("physics") ;
  sample.defineType("control") ;

  // Construct combined dataset in (x,sample)
  // 将数据导入到RooDataSet里面,其中有目录Index(sample),还有变量x,并将目录中的名字与数据集关联起来
  RooDataSet combData("combData","combined data",x,Index(sample),Import("physics",*data),Import("control",*data_ctl)) ;

  // 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 )
  // -----------------------------------------------------------------------------------
  // Construct a simultaneous pdf using category sample as index
  // 这里我们将之前定义的目录送入RooSimultaneous类中去,然后
  RooSimultaneous simPdf("simPdf","simultaneous pdf",sample) ;

  // 这里将之前定义的分布模型放到RooSimultaneous类的sample里面去,
  // Associate model with the physics state and model_ctl with the control state
  simPdf.addPdf(model,"physics") ;
  simPdf.addPdf(model_ctl,"control") ;

  // 这里将数据进行拟合,拟合的目标物理分布的模型就是前面的
  // P e r f o r m   a   s i m u l t a n e o u s   f i t
  // ---------------------------------------------------
  // Perform simultaneous fit of model to data and model_ctl to data_ctl
  simPdf.fitTo(combData) ;
  // 到此为止,就完成了拟合的计算过程,其实已经结束了,只是不直观要画图
  // 可以看出,RooSimultaneous就是一个类中包含了多个模型,这些多个模型的公共参数用来描绘共同的一个物理分布。
  // 这种多个模型,可以来自于不同衰变道,也可以来自于不同的样本如例子中随机产生的物理分布和控制样本。

  // 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
  // ----------------------------------------------------------------
  // Make a frame for the physics sample
  // 画图的RooPlot,用来载入画图的内容,在这里定义了分bin
  RooPlot* frame1 = x.frame(Bins(30),Title("Physics sample")) ;

  // Plot all data tagged as physics sample
  // 将数据画在上面定义的frame1中,并且通过Cut条件只看physics目录的内容
  combData.plotOn(frame1,Cut("sample==sample::physics")) ;

  // Plot "physics" slice of simultaneous pdf.
  // NBL You _must_ project the sample index category with data using ProjWData
  // as a RooSimultaneous makes no prediction on the shape in the index category
  // and can thus not be integrated
  simPdf.plotOn(frame1,Slice(sample,"physics"),ProjWData(sample,combData)) ;
  // 这里可以看到是又画了一次,这次是单独画的本底用的kDashed的线,可以按后面给的方法搜索类型继续看
  simPdf.plotOn(frame1,Slice(sample,"physics"),Components("px"),ProjWData(sample,combData),LineStyle(kDashed)) ;

  // 定义第二个图的画图的frame
  // The same plot for the control sample slice
  RooPlot* frame2 = x.frame(Bins(30),Title("Control sample")) ;
  combData.plotOn(frame2,Cut("sample==sample::control")) ;
  simPdf.plotOn(frame2,Slice(sample,"control"),ProjWData(sample,combData)) ;
  simPdf.plotOn(frame2,Slice(sample,"control"),Components("px_ctl"),ProjWData(sample,combData),LineStyle(kDashed)) ;

  // 定义好画布,将前面的frame调用Draw()成员函数,就画图了
  TCanvas* c = new TCanvas("rf501_simultaneouspdf","rf403_simultaneouspdf",800,400);
  c->Divide(2) ;
  c->cd(1) ; gPad->SetLeftMargin(0.15) ; frame1->GetYaxis()->SetTitleOffset(1.4) ; frame1->Draw() ;
  c->cd(2) ; gPad->SetLeftMargin(0.15) ; frame2->GetYaxis()->SetTitleOffset(1.4) ; frame2->Draw() ;
}

实际上,我们看完上面的内容就可以了解到其实我们什么都没有做!我们做的仅仅是:

头文件
函数{
  定义预期的多个模型model(RooAddPdf)
  数据data(DataSet),按照样本目录加载
  调用联合拟合类(RooSimultaneous).fitTo(data)
  画图,信号本底
}

其他的当然也是类似,是计算机帮我们完成了这些操作,当然如涉及到其背后是如何工作的则要继续了解simultaneous背后的数学是如何工作的,比如具体的误差它是如何计算的?整体上参数是如何传递,计算是如何迭代的?编程?不存在的,可以放后去了解。这也学就是我们称之为脚本的原因吧,之所以使用脚本文件,就是相应的程序包提供了简易的借口,可以对应用的用户功能进行自定义。那么有没有办法不写这么多字呢好头大啊?有的,有的人用Python来自动写脚本,有的人直接用PyROOT。

使用Vim

如何使用vim?

这个实际上吧vim的手册简单看一下,然后了解一下计算机发展过程,“文本”这个概念是如何发展的,就可以更好地使用vim,我自己在两年用了杨振伟老师推荐的一个简单的vimrc文件易水博客,还有重新编译了vim用来安装了vim-snippets后,再没有怎么搞过了,我想snippets应该已经功能很强大了。 想起来第一次用vim其实就和我小时候玩游戏一样,上来就键盘挨个按一遍,软件也是,把菜单打开所有按钮按一遍,这就是我从来看任何计算机相关东西的所谓“窍门”。

继续试着学ROOT有什么方便的方法?

如果有机会找人教,可以问一下别人,或者有人给脚本,就可以学的很快,因为大家公用的东西都是比较成熟的。在简单用了之后,还要继续学懂一些,或者看到别的ROOT脚本文件有时候发现看不懂,甚至受到了惊吓,下面的一点内容我们试着想一想办法。

首先想到的就是,看一看例子tutorial,ROOT官方带有许多的例子,只要稍加更改就可以学到很多东西,上面篇幅中做联合拟合的内容就是依据这样的方式来做的,下面还有。

一个是搜索使用的类的定义和继承关系,在网站这里ClassIndex,比如可以试着搜索ProjWData这个成员,也可以在ClassList用搜索框或者是按类找名字ctrl+f,这样搜索可以找到ROOT编程的详细内容。 可以看到,其定义有两个如下:

第一个定义可以看出,其可有两个变量,一是数据集RooDataSet(RooAbsData三个成员中的一个),从下面的继承图可以清楚地了解到程序的调用关系,第二个可能用于分bin的问题,

RooFit::ProjWData  ( const RooAbsData & projData,
Bool_t  binData = kFALSE
)

从这里上下两个定义可以看出,ROOT软件包整体内部命名是十分规范的,我们在使用的时候也要尽量把变量函数名写的清楚易懂有一定的规范。如BES合作组的BOSS分析环境,可以查看BOSS定义的基本的类和继承图BOSS Class Reference Documentation,这里可以清楚看到很多分析中使用的内容。当然,如果有较为统一的规范是不容易的,因为在传统研究中各种不同群体都有可能产生各自内部的“独特语言系统”。

第二个定义可以看出,其可有三个变量,一是从数据集里取数据的方法(RooCategory是其中一个),第二个见上,第三个可能用于分bin的问题

RooFit::ProjWData  ( const RooArgSet & projSet,
const RooAbsData & projData,
Bool_t binData = kFALSE
)

可以看出我们这里使用了两次,一次第一个定义,一次第二个定义,并且没有用bin变量。 这样我们通过查看这些的定义关系,有助于梳理清楚他们的调用关系,用起来就不会发怵。 所谓计算机程序就是一大堆的相互调用,从ifelse编译一直调用到逻辑门基本不会出错,如果出现问题,就查看定义看哪里调用跑丢了,如果能够解决好把调用的类型放对了,一般程序都能正常运行,就像喂小猫一样。 如果想要训练一下编程时候的脑子,其实有更直观一点的编程比如说用G语言LabVIEW,你可以看数据流,并且简单的编程结构也比较直观,实际应用也比较广泛,比如BESIII值班室的用户界面系统,可以看到很多都是使用LabVIEW编的。另外甚至可以通过一些软件也能训练,如Zachronic公司旗下的软件都很硬核,还有辅助学习各种编程语言的网站酷代码非常有趣,很全各种计算机语言都能找到。

可以在网上ROOT Forum看别人遇到了什么问题,如果身边没人也可以试着问,还可以看看同行都在关心什么问题? 这里提供一个问问题模板参考:

Hi,
I want to do sth.
I 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.
I need to to do similar processing but for a different "sth.? dataset?".
Could some body guide me to the Documentation or Tutorials for doing this?

Thanks
"yourname"

Sample: "your sample idea"
Code:"your sample code"

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

~~继续思考~~

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