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) ;

//RooAddPdfmodel("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用作峰位的初始化值很方便，且很准确