专栏首页人芳觅ROOT-数据读取-直方图-Roofit拟合基本流程-(入门实用)

ROOT-数据读取-直方图-Roofit拟合基本流程-(入门实用)

笔者最近在测核素能谱,测出的能谱需要分析,比如计算某全能峰的分辨率。用到的数据处理分析工具是ROOT(cern),整个能谱读取分析的流程可给各位看官当入门或干货材料使用。不过,ROOT大神就必看本文了,至少节约2分钟的时间,日后要是有新鲜的ROOT技巧会另作分享。

本文ROOT分析的实验背景和数据处理目的如下图1所示:

图1 实验背景和ROOT数据处理的目的

目的1:得到扣除本底后的能谱并绘制出图

基本过程为:1)*.Spe文件数据读取;2)定义一个ROOT文件用于存储能谱;3)数据读取到直方图,并作本底扣除;4)写入ROOT文件;5)画图。

#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooBreitWigner.h"
#include "RooPolynomial.h"
#include "RooAddPdf.h"
#include "RooFitResult.h"
#include "RooRealIntegral.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "TFile.h"
#include "TTree.h"
#include "TH2D.h"
#include "TH1D.h"
#include "THStack.h"
#include "TROOT.h"
#include "TLatex.h"
#include "TF2.h"
#include "TF1.h"
#include "iostream"
#include "fstream"
#include "RooCBShape.h"
#include "TMath.h"
#include "RooFit.h"
#include "RooMath.h"
 using namespace RooFit;
 using namespace std;
void templateread(){
const int varnum=2;
const char *newlabelname[varnum]={"CsI-B","CsI-Cs137"};
string filename[varnum]={"20190401-csi-bg-1725s.Spe","cs137-csi-1725s.Spe"};
double measuringtime[varnum]={1725,1725};
       const int channelnum=24+2048;
       string tempdata;
       ifstream fin[varnum];
       TH1F *spec[varnum];
 
int font=22;
TFile *outputFile=new
TFile("CsI-1cm1cm1cm-new.root","RECREATE");
for (int i1=0;i1<2;i1++)
{      fin[i1].open(filename[i1],ios_base::in);
        spec[i1] = new TH1F(newlabelname[i1],newlabelname[i1],2048,0,2048);
       for (int i=0;i<channelnum;i++)
       {
                fin[i1]>>tempdata;
                if(i>24)
 
                {
                       spec[i1]->SetBinContent(i-24,stod(tempdata));
                }
 
}
 
spec[i1]->SetFillColor(kRed);
spec[i1]->SetYTitle("Counts");
spec[i1]->SetXTitle("Channel");
spec[i1]->SetLineColor(kRed);
spec[i1]->SetLabelSize(0.04,"X");
spec[i1]->SetLabelSize(0.04,"Y");
spec[i1]->SetLabelFont(font,"X");
spec[i1]->SetLabelFont(font,"Y");
spec[i1]->SetLabelOffset(0.01,"X");
spec[i1]->SetLabelOffset(0.01,"Y");
spec[i1]->SetTitleOffset(1.0,"X");
spec[i1]->SetTitleOffset(1.2,"Y");
spec[i1]->SetTitleOffset(1.2,"Z");
spec[i1]->SetTitleSize(0.05,"X");
spec[i1]->SetTitleSize(0.05,"Y");
spec[i1]->SetTitleSize(0.05,"Z");
spec[i1]->SetTitleFont(font,"X");
spec[i1]->SetTitleFont(font,"Y");
spec[i1]->SetTitleFont(font,"Z");
 
}
 
gStyle->SetOptStat("");
 
spec[1]->Add(spec[0],-(measuringtime[1]/measuringtime[0]));
for (int j=0;j<2;j++)
{
   for (int i=0;i<2048;i++)
    {
 
       if(spec[j]->GetBinContent(i+1)<0)    {spec[j]->SetBinContent(i+1,0);}
 
    }
}
for (int i=0;i<20;i++)
{
spec[1]->SetBinContent(i+1,0);
}
outputFile->Write();
//spec[0]->Write();
//spec[1]->Write();
//outputFile->Close();
TCanvas*myc=new
TCanvas("myc","myc",650,450);
myc->SetLeftMargin(0.12);
myc->SetRightMargin(0.08);
myc->SetTopMargin(0.08);
myc->SetBottomMargin(0.12);
spec[1]->Draw();
//myc->Print(Form("CsI-NoBackground-sec1725.png"));
}

图2 绘制能谱

目的2:对能谱全能峰进行拟合

基本过程为:1)已经保存的ROOT文件中读取直方图;2)定义拟合函数的参数区间;3)选择感兴趣的几个函数用于全能峰拟合;4)绘制拟合结果。

#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "TAxis.h"
//#include <stdlib.h>
using namespace RooFit ;
using namespace std ;
 
void goodfit()
{
//TFile f1("20190501bgona22.root");
//TH1D* hist =(TH1D*)f1.Get("BGO-Na22");
TFile* inputFile = TFile::Open("20190501bgona22.root");
//TFile* inputFile = TFile::Open("CsI-1cm1cm1cm-new.root");
//TFile* inputFile = TFile::Open("GAGG-6mm6mm6mm-cs137-exp.root");
TH1F *spec[4];
//spec[0]=new TH1F();
spec[0]=(TH1F*)inputFile->Get("BGO-Na22");
//spec[0]=(TH1F*)inputFile->Get("Cs137-TC");
//spec[0]=(TH1F*)inputFile->Get("CsI-Cs137");
//spec[0]=(TH1F*)inputFile->Get("BGO-Cs137");
// S e t u p   m o d e l
// ---------------------
//double myscale=1.0/h1->Integral();
//h1->Scale(myscale);//normalize the hist
int meanv=450;
int rangexmin=390;int rangexmax=580;
int xmin=rangexmin;int xmax=rangexmax;
RooRealVar x("x","x",xmin,xmax) ;
RooRealVar mean("mean","mean",meanv,rangexmin,rangexmax) ;
RooRealVar sigma("sigma","sigma",20,5.,500) ;
RooGaussian sig("sig","signal p.d.f.",x,mean,sigma) ;

RooRealVar argpar1("argpar1","argus shape parameter",20,0.,40.) ;
RooRealVar argpar2("argpar2","argus shape parameter",20,0.,40.) ;
RooArgusBG argus("argus","Argus PDF",x,argpar1,argpar2) ;
RooRealVar a0("a0","a0",-0.0001,-1.,0.1);
// RooRealVar a1("a1","", 0.5, -1, 100);
RooExponential bkg("bkg","background p.d.f.", x,a0);
RooRealVar c0("c0","coefficient #0", -1.0,-300.,10.) ;
RooRealVar c1("c1","coefficient #1", 0.1,-100.,1.) ;
RooRealVar c2("c2","coefficient #2",-0.1,-100.,2.) ;
RooChebychev compton("bkg","background p.d.f.",x,RooArgList(c0,c1,c2)) ;
RooRealVar mean_bkg("mean_bkg","mean",rangexmin/2+rangexmax/2,rangexmin,rangexmax);
RooRealVar sigma_bkg("sigma_bkg","sigma",30,10.,60.);
RooGaussian bkg_peak("bkg_peak","peaking bkg p.d.f.",x,mean_bkg,sigma_bkg) ;
RooRealVar fpeakbkg("fpeakbkg","peaking background fraction",0.5,0.4,1.) ;
RooAddPdf totalbkg("totalbkg","compton+bkg_peak",RooArgList(bkg_peak,compton),fpeakbkg);
 
RooRealVar fsig("fsig","signal fraction",0.9,0.5,1.) ;
//RooAddPdf
model("model","sig+(compton+bkg_peak)",RooArgList(sig,totalbkg),fsig);
RooAddPdf model("model","sig+bkg-e",RooArgList(sig,bkg),fsig) ;
//RooAddPdf
model("model","sig+compton",RooArgList(sig,compton),fsig) ;
//RooAddPdf
model("model","sig+compton",RooArgList(sig),fsig) ;
//RooAddPdf
model("model","sig+compton",RooArgList(sig,argus),fsig) ;
 
//RooPlot* frame = x.frame(Title("CsI Cs-137 662keV Photopeak Fitting;Channel"));
//RooPlot* frame = x.frame(Title("BGO Cs-137 662keV Photopeak Fitting;Channel"));
RooPlot* frame = x.frame(Title("BGO Na-22 511keV Photopeak Fitting;Channel"));
//RooPlot* frame = x.frame(Title("GAGG Ba-133 81keV Photopeak Fitting"));
RooAbsData* data = new RooDataHist("data","data",x,spec[0]);
data->plotOn(frame);
//data->plotOn(xframe,Binning(1000));
//model.fitTo(*data,Save(),Extended());
model.fitTo(*data,Save(),Extended(),Range(rangexmin,rangexmax));
//model.fitTo(*data,Save(),Extended(),Range(rangexmin,rangexmax));
 
model.plotOn(frame,LineColor(kGreen)) ;
//model.plotOn(frame,Components(argus),LineColor(kBlue),LineStyle(kDashed)) ;
model.plotOn(frame, Components(compton),LineColor(kBlue),LineStyle(kDashed)) ;
model.plotOn(frame, Components(sig),LineColor(kRed),LineStyle(kDashed)) ;
//model.plotOn(frame, Components(bkg_peak),LineColor(kViolet),LineStyle(kDashed)) ;
//model.plotOn(frame, Components(totalbkg),LineColor(kBlue),LineStyle(kDashed)) ;
//model.plotOn(frame, Components(bkg),LineColor(kBlue),LineStyle(kDashed)) ;
 
//model.paramOn(frame, Format("NEU"), Layout(0.15,0.50,0.9));
model.paramOn(frame,Parameters(RooArgSet(sigma,mean)), Format("NEU"), Layout(0.55,0.9,0.9));
//model.paramOn(frame,Parameters(RooArgSet(sigma,mean,c0,c1,c2)), Format("NEU"),Layout(0.55,0.9,0.9));
//model.paramOn(frame,Parameters(RooArgSet(sigma,mean,fsig,a0)), Format("NEU"), Layout(0.55,0.9,0.9));
TCanvas* c = new TCanvas("c2","Fitting test",1200,800);
frame->Draw() ;
//c->Print(Form("Na22-511keVfit0903.png"));
//c->Print(Form("Ba133-81keVfit0903.png"));
}

图3 全能峰拟合

到这里,两个目的均已达成,Roofit其实算是种很偷懒的拟合,未来的教程将探讨普适的Fit以及TSpectrum的机智用法。

图4 TSpectrum用作峰位的初始化值很方便,且很准确

喜欢的话,分享一下吧~^o^~

本文分享自微信公众号 - 人芳觅(lililiangzuozuo),作者:梁佐佐

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2019-05-13

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • Geant4--一次编译,运行多个Run,极大提升模拟效率

    应唐光毅博士/后之约,对于Geant4模拟,我们看是否能解决这么一个问题:我现在想模拟探测器不同角度下的响应,每次模拟需要/run/beamOn 100, 可是...

    梁佐佐
  • 回答某一个同学Geant4问题;建立了一个Geant4群聊,欢迎大家来踩

    要想熟练得使用Geant4,C++的基本知识是要掌握的,其学习参考教材可任意选,刚入门的同学推荐比较薄的谭浩强版本《C++面向对象程序设计》,网上资源也丰富,比...

    梁佐佐
  • 3D建模(Blender)+实例演示透明材质和5毛光晕特效

    不知道为什么,明明只想安安静静的科研(搬铅砖,一块22斤)。却发现自己,不会P图是不行的,毕竟你总想在自己的PPT里面突出点什么;不会视频剪辑、不会Gif动图,...

    梁佐佐
  • 如何在保留装箱对象的前提下修改值

    有人问如何在保留装箱对象的前提下修改值? 场景: object obj = 100; Console.WriteLine("original object va...

    葡萄城控件
  • 给Win10瘦身

    Windows 10自带的应用程序,其中包括时钟、日历、地图、天气、相册等多款应用既占硬盘空间又没什么用。那个新的照片查看器界面又丑又难操作,影响我欣赏小姐姐的...

    fem178
  • D3D还是GDI? Windows平台播放RTSP或RTMP渲染模式比较

    好多开发者纠结播放端绘制,是D3D还是GDI,先说结论,Windows平台播放渲染这块,支持D3D的前提下,优先D3D,如果检测到不支持D3D,数据回调上来,G...

    daniulive
  • 《黑盾》首部VR宣传片曝光震撼视角斩杀好莱坞大片

    VRPinea
  • [菜鸟SpringCloud入门]第二章:创建服务提供者并在Eureka进行注册

    欢迎来到菜鸟SpringCloud入门实战系列(SpringCloudForNoob),该系列通过层层递进的实战视角,来一步步学习和理解SpringCloud。

    Rude3Knife的公众号
  • hdu1095

    @坤的
  • CSDN 阅读 源代码 ListView宽度过小获取博客

    我们写软件,一般是用用一个页面来做导航,这个页面就是一个Frame,然后包含各种导航,所以这个页面会一直存在我们的内存。

    林德熙

扫码关注云+社区

领取腾讯云代金券