Table of Contents

• 1  一、ROOT软件操作模式介绍
  • 1.1  1.1 终端模式常用命令
  • 1.2  1.2 对话框模式
  • 1.3  1.3 脚本模式
• 2  二、脚本模式相关功能代码介绍
  • 2.1  2.1 创建带有tree结构的root文件大致流程
  • 2.2  2.2 读取ROOT的tree数据，逐事件分析流程
  • 2.3  2.3 头文件
  • 2.4  2.4 C++/ROOT常用字符串输入输出
    • 2.4.1  2.4.1 代码内输入输出
    • 2.4.2  2.4.2 将字符存储成文本
    • 2.4.3  2.4.3 主函数参数argc与argv说明
  • 2.5  2.5 MakeClass方法
    • 2.5.1  2.5.1 调用MakeClass方法
    • 2.5.2  2.5.2 MakeClass with TChain
  • 2.6  2.6 探测器刻度（扣本地、寻峰、排序对应、拟合、提取参数）
    • 2.6.1  1 扣本底--TSpectrum (选择使用)
    • 2.6.2  2 寻峰--TSpectrum
    • 2.6.3  3 峰位排序-map
  • 2.7  2.7 在root文件中创建/读取子文件夹
    • 2.7.1  2.7.1 创建
    • 2.7.2  2.7.2 读取
  • 2.8  2.8 TCut的使用
    • 2.8.1  2.8.1 ROOT终端调用
    • 2.8.2  2.8.2 TCut存成脚本，在终端调用
    • 2.8.3  2.8.3 TCut存成脚本，在脚本中调用，判断是否在区域内；
• 3  三、画图操作命令（画图、拟合、坐标修改、颜色更换等等）
  • 3.1  3.1 脚本画图
  • 3.2  3.2 得到histogram或TGraph指针（两种方法注意区分）
    • 3.2.1  3.2.1 得到Histogram指针
    • 3.2.2  3.2.2 得到TGraph指针
    • 3.2.3  3.2.3 Histogram转TGraph，以及命名
  • 3.3  3.3 叠图
  • 3.4  3.4 一个画布画多个图
    • 3.4.1  方法1
    • 3.4.2  方法2
  • 3.5  3.5 填加图例
  • 3.6  3.6 histogram拟合与结果提取
  • 3.7  3.7 脚本存储pdf文件
  • 3.8  3.8 其他功能
  • 3.9  3.9 颜色以及功能项
    • 3.9.1  3.9.1 Colors
    • 3.9.2  3.9.2 width and style of Line
    • 3.9.3  3.9.3 width and style of Marker
    • 3.9.4  3.9.4 option of tree draw and graph
    • 3.9.5  3.9.5 option of TGraph fit
    • 3.9.6  3.9.6 option of TH1 fit
• 4  四、功能代码
  • 4.1  4.1 利用TMath::Sort()对能量信号排序
    • 4.1.1  4.1.1 对于Hit结构数据
    • 4.1.2  4.1.2 对于非Hit结构数据
  • 4.2  4.2 同一粒子在硅条正背面响应的挑选
  • 4.3  4.3 由硅条正背面响应条(xstrip,ystrip)转换成($\theta, \phi$)
  • 4.4  4.4 同一望远镜两层硅条的PID匹配
  • 4.5  4.5 激发能谱+效率曲线(放大)
  • 4.6  4.6 激发能谱本底拟合(TH1D*)->ShowBackground()

一、ROOT软件操作模式介绍

• 本内容纯为个人爱好总结，相关内容以root官网为准https://root.cern/doc/master/index.html
• ROOTPrimer and ROOTUserGuide下载(提取码：8zm8)

1. 终端模式（进入root软件的terminal通过键入有关命令简单的查看、检验数据）
2. 对话框模式（进入root软件后，键入TBrowser b）
3. 脚本模式（通过编写脚本文件，实现数据的批量分析、处理、画图、存图等功能）

注意

• 脚本模式指令与终端模式命令基本相似，但脚本模式某写功能比终端模式更严格，需要作些准备工作！
  • 比如终端模式tree->Draw()即可实现画图；脚本模式测需要创建画布、填图、画图；
• 终端模式与对话框模式下，对tree->Draw()的操作可通过鼠标实现；
• 脚本模式下，对tree->Draw()的操作仅能键入对应的命令才能实现；

1.1 终端模式常用命令

(对于tree->Draw()的具体操作单独介绍)

• 打开root文件

root -l path+name.root

• 查看当前root中结构信息：root文件名字、tree结构名字、存储histgram-1D/2D/3D等信息

.ls

• 查看tree结构信息，即tree中Branch具体内容,为后面画图作准备
  • 此处tree为.ls执行之后，TTree之后的tree的名字，此处的tree就叫tree；

tree->Print();

• 读取tree中第event number个事件的信息；

tree->Show([event number])

• 浏览tree中第event number个事件的信息；

tree->Scan("Branchname1:Branchname2:...:Branchnamen","","",numbers of events, start events);

• 画tree结构中Branch的histogram
  • 如tree->Draw("branch_y:branch_x>>alias(xbin,minx,maxx,ybin,miny,maxy)","condition,cut etc.","colz etc.");
  • 第一个双引号：内容+座标轴范围与bin，alias为该条件别名；
  • 第二个双引号：画图的条件、符合、CUT等；
  • 第三个双引号：option，如colz等；

tree->Draw("ctof:td-tu>>(2000,-20,50,1000,0,100)","PID==0","colz");
tree->Draw("td-tu>>htx(500,-20,50)","","");

• 画root中以画布形式存储（非tree结构）histogram图；

hctof->Draw();//hctof可从.ls命令查看

• excute outer command in terminal of root

  • ![shell command]
  • eg:!ls -lh treeADC.*

• 退出root软件

  • .q

1.2 对话框模式

(进入后有关拟合、cut、坐标轴、曲线等功能可通过鼠标实现)

• 进入对话框模式

TBrowser b

1.3 脚本模式

利用C++语言，借由ROOT库函数，自己编写代码（.c/.cpp/.h等等）用于数据分析处理以及批量画图、存图等；

• 编写分析代码并命名：如filename.c;
• 在终端键入如下命令：

root -l filename.c

二、脚本模式相关功能代码介绍

2.1 创建带有tree结构的root文件大致流程

1. 主函数体(文件名要与主函数名称相同)；
2. 常量声明，一般包括物理常量/探测器常量；
3. tree结构中Branch所需的变量声明；
4. 定义新的root文件，声明新的tree；
5. 将变量地址添加到tree中，其tree中Branch与变量地址关联；还可以定义hisrogram；
6. 循环体，计算各种变量，并填入数据；
7. 将数据写入histogram/tree，存储并关系root文件；

详细代码见李老师数据分析课homework1.1中的tree.cc

代码示例

////主函数声明
void/Int_t main(){}

////变量声明
Int_t a;
Double_t b;
Long64_t c;
......

////创建名为tree.root的文件，指针为*opf,与opf->Close()匹配使用，告知程序.root文件处理完成/关闭；
TFile * opf = new TFile("tree.root","recreate");
////创建新tree，名称为tree，其指针为*opt;
TTree * opt = new TTree("tree","tree structure");

////为opt指针创建Branch("变量名称",变量指针,"变量名/变量数据格式")
////将之前声明的变量指针与Branch链接一起，注意此处Branch相关变量要与之前声明的变量对应；
opt->Branch("x",&x,"x/D");
////按照opt->Branch(" ",  ," ")(一般为多个Branch)所定义的参数向opt指针填入全部数据；
opt->Fill();//该命令一般放入循环体中
////创建histogram,TH1D、TH2D、TH3D、TH1F、TH2F、TH3F等
TH1D * hctof = new TH1D("hctof","neutron time of flight",1000,0,100);
////将数据填入刚刚声明的histogram中，Fill(x,y,z)    
hctof->Fill(ctof);//该命令一般放入循环体中

////将声明的histogram写入root文件；
hctof->Write();
////将tree中Branch有关数据写入root文件；
opt->Write();
////关闭root文件
opf->Close();

2.2 读取ROOT的tree数据，逐事件分析流程

！！！注意输入tree与输出tree的差异，特别是Makefile编译方法时！！！

1. 主函数；
2. 打开root文件，得到tree的TTree指针；
3. 声明tree结构中Branch所需的变量；
4. 将tree中Branch的数据的指针指向声明的变量；
5. 将数据写入新的ROOT文件（参照《创建带有tree结构的root文件》）；
6. 得到tree的事件总数，由循环体遍历每个事件；
7. 将第jentry个事件数据填入Branch中指向的变量，填入对应的fistogram与tree；
8. 画图/关闭tree，写入新的tree，关闭新的root，关闭老的root；

详细代码见李老师数据分析课homework1.1中的homework1.1-sum

代码示例

//打开名为tree.root的文件（只读），并指定指针为ipf，脚本最后注意关闭文件ipf->Close();
TFile * ipf = new TFile("tree.root", "read");
//检验tree.root文件是否存在；
if (ipf->IsZombie())
{
    cout << "Error opening file" << endl;
    exit(-1);
}
//为了避免多个root数据扰乱，加入此代码;
ipf->cd();
//得到名字为tree的TTree指针，tree的名字可由.ls得到
TTree * tree = (TTree*)ipf->Get("tree");
//声明tree的Branch变量（可由tree->Print()得到或通过MakeClass命令生成相关文件）
Double_t x;
Int_t pid;
//将变量指向tree中Branch的地址(注意此处与创建Branch格式不同，没有数据类型）
tree->SetBranchAddress("x",&x);
tree->SetBranchAddress("pid",&pid);
//---将新数据写入新的ROOT文件（参考创建带有tree结构的root文件大致流程）
//---新变量声明；
Double_t x_cali;
//---创建新root、新tree；
TFile * opf = new TFile("tree.root","recreate");
TTree * opt = new TTree("tree","tree structure");
//---新变量地址指向新的tree的Branch；
opt->Branch("x_cali",&x_cali,"x_cali/D");

//逐事件读取tree的branch数据
//得到tree的事件总数，利用for循环遍历事件
//if (fChain == 0) return;
Long64_t nentries = tree->GetEntries();//！！！此处tree为输入文件的tree！！！
Long64_t nbytes = 0, nb = 0;
for(Long64_t jentry=0;jentry<nentries;jentry++)
{
//循环体中提取第n个事件的信息(在此之前已经声明了该root文件中的tree的Branch与变量地址的对应关系)；
    tree->GetEntry(jentry);//！！！此处tree为输入文件的tree！！！
    Long64_t ientry = LoadTree(jentry);
    if (ientry < 0) break;
    nb = tree->GetEntry(jentry);    nbytes += nb;
//---新的变量 = 旧数据处理后结果（参考创建带有tree结构的root文件大致流程）
    x_cali=f(x);
//---将新的结果填入tree中；（参考创建带有tree结构的root文件大致流程）
    opt->Fill();////！！！此处tree为输出文件的tree！！！

//为了检查数据处理进程，可添加
    if(jentry%100000==0) cout<<"process"<<jentry<<"of"<<nentries<<endl;
}

//关闭读入的root文件
ipf->Close();
//---将tree中Branch有关数据写入root文件；
opt->Write();
//---关闭root文件
opf->Close();

2.3 头文件

自己编写可编译执行文件时

• <>:为ROOT、C++库函数，
• ""：为自己编写的.h,.c等函数
• C++(相关函数头文件请查看：http://www.cplusplus.com/)

#include <iostream> //涉及cout，cin
#include <sstream>  //sprintf、string
using namespace std;    //std标准输入输出

//ROOT(相关函数头文件请查看：https://root.cern/doc/master/index.html)
#include <TROOT.h>  //ROOT标准
#include <TChain.h> //Chain：多个root首位连结（数据结构要求一样）
#include <TFile.h>  //TFile：打开、创建.root文件
#include <TTree.h>  //TTree *tree:创建tree结构；
#include <TH2.h>    //TH2I,TH2D,TH2F:创建二维histogram；
#include <TStyle.h> //gStyle->SetPalette()
#include <TCanvas.h>    //TCanvas * c1 = new TCanvas(""):创建画布；
#include <TF1.h>    //pol1、gaus，TF1 * f1:调用函数以及创建函数；
#include <TFitResult.h> //拟合、调整参数、提取拟合：参数
#include <TString.h>    //sstream的root化；
#include <TGraph.h> //TGraph * g1:创建graph，注意区graph vs histogram；

条件编译 #ifndef _HEAD_H...中下划线的理解

• 假设你的头文件名为norm.h，根据习惯，我们声明一个宏NORM_H对应这个头文件。在头文件中开始的地方和结尾的地方加上#ifndef NORM_H以及#define NORM_H实现对HEAD_H的声明和判断，头文件norm.h如下：

#ifndef   NORM_H 
#define   NORM_H 

//填写nrom.h具体内容

#endif

• __FILE__,__LINE__ 都是与定义的宏，使用“_”和“__ ”开始的函数一般都是专用的函数，一般都是于特定系统相关的，如果要想有更好的移植性，应该避免使用。

2.4 C++/ROOT常用字符串输入输出

2.4.1 代码内输入输出

//c++
sprintf(charname,"c1data%04d.root",i);
string path = "/home/wangdongxi/ana_xi/caliprepare/c2forcaliroot/";

//root
//声明
TString name;
//赋值
name = "内容"; 
name.Form("内容1%s%d内容2", 变量1, 变量2); // %s/%d为变量格式, 具体可见格式化输出
//调用
name.Data()

//声明+赋值
TString name1 = TString::Format("内容1%s%d内容2", 变量1, 变量2);  //%s/%d为变量格式, 具体可见格式化输出;

2.4.2 将字符存储成文本

ofstream writer0d;
writer0d<<fixed<<setprecision(10);//设置输出精度；
TString outputr0d =TString::Format("./result_sum/r0d_result.txt");
writer0d.open(outputr0d.Data(),ios::out);
writer0d<<setw(15)<<"runnumber"<<" "<<setw(15)<<"sigma"<<endl;
for(auto im = r0dsigma.begin(); im != r0dsigma.end(); im++)
{
     writer0d<<setw(15)<<im->first<<" "<<setw(15)<<im->second<<endl;
}
writer0d<<setw(15)<<"runnumber"<<" "<<setw(15)<<"mean"<<endl;
for(auto im = r0dmean.begin(); im != r0dmean.end(); im++)
{
    writer0d<<setw(15)<<im->first<<" "<<setw(15)<<im->second<<endl;
}
writer0d.close();

2.4.3 主函数参数argc与argv说明

给主函数传递参数时，可用下面方法：

int main(int argc, char** argv)
{
//1 argc：提供给主函数的参数个数；
//2 argv：参数的字符串数组指针，其中argv[0]指向程序运行的全路径名，argv[1]指向程序名的第一个字符串，argv[2]指向程序名的第二个字符串，……；
//3 argv是字符串类型，如需整型数据可用atoi(argv[1])转换数据类型；
}

2.5 MakeClass方法

作用

• 自动生成该root文件的数据结构，并生成循环执行脚本（Loop()函数）；
• 直接添加想要的内容，并在Loop()中添加分析程序，无需从头编写（如2.2节）；
• 对于编写可编译执行程序，采用继承的方法链接该文件，免去了键入输入文件tree结构的内容；

详细方法见李老师数据处理课homework2.1中的exercise2.4以及作业分析代码

2.5.1 调用MakeClass方法

tree->MakeClass("ppac");//输出ppac.h和ppac.C文件
//修改ppac.C中Loop部分后
gROOT->ProcessLine(".L ppac.C");//ROOT命令行.L ppac.C
ppac t;//实体化
t.GetEntry(13);//运行ppac类中GetEntry()函数
t.Show();//运行ppac类中Show()函数
t.Loop();//运行ppac类中Loop()函数,逐个事件分析数据；

2.5.2 MakeClass with TChain

TChain* chain=new TChain("tree");//chain即位单个root文件中的tree
chain->Add("run0001.root");
chain->Add("run0002.root");
chain->Add("run0003.root");
chain->Add("run0004.root");
chain->Add("run0005.root");
chain->MakeClass("ppac_chain");//输出chain脚本文件

2.6 探测器刻度（扣本地、寻峰、排序对应、拟合、提取参数）

具体代码请看李老师数据处理课homework3.1中的拟合与刻度

1 扣本底--TSpectrum (选择使用)

//用法：Int_t TSpectrum::Background(const TH1 * hin,   //input 1-d histogram
//                                                  Int_t niter = 2,        //光滑度
//                                              Option_t * option = “”) //功能项

TSpectrum *s=new TSpectrum(500);                //实例化，500：maximum number of peaks
TH1F *hb=(TH1F*)s->Background(h0,30,"same");        //拟合h0本底，并传递给hb，30-本底光滑程度, 用户进行调节
h0->Clone("hh");        //克隆h0，命名为hh
hh->Add(hh,hb,1,-1);        //将hh与hb叠加，TH1F * h0 = Add(TH1F * h1, TH1F * h2, c1, c2)
//                                       即h0 = h1 * c1 + h2 * c2

2 寻峰--TSpectrum

//用法：Int_t TSpectrum::Search(const TH1 * hin,       //input 1-d histogram
//                  Double_t sigma = 2,     //光滑度
//                  Option_t * option = “”,     //功能项
//                  Double_t threshold = 0.05)  //count>ther * Max of Counts of peak
Int_t nfound=s->Search(hh,2,"",0.05);   //寻峰的结果（个数）传递给nfound；
Double_t *xpeaks,*ypeaks;           //指针
xpeaks=s->GetPositionX();           //寻峰的结果（Mean）传递给xpeak；
ypeaks=s->GetPositionY();           //寻峰的结果（Counts）传递给ypeak；
//A pointer to the polymarker object can be retrieved later via:
TList *functions = hin->GetListOfFunctions();
TPolyMarker *pm = (TPolyMarker*)functions->FindObject("TPolyMarker");

3 峰位排序-map

//  ***map是STL（中文标准模板库）的一个关联容器。
//  ***map提供一对一的数据处理，key-value键值对，其类型可以自己定义，第一个称为关键字，第二个为关键字的值。
//  ***multimap可以允许多个key值是相同的，而map的key值不允许重复。
//  ***map,multimap, 按照key值的大小自动排序，key值有小到大。
//  ***需要使用头文件<map><alforithm>

//定义map变量
map<Int_t,Double_t> m1;//第一个是键（key）的类型，第二个是值（value）的类型
multimap<Int_t,Double_t>mul1;
//赋值
TRandom3 *g=new TRandom3(0);
for(int i=0;i<nfound;i++)
{
    int key=xpeaks[i];
    double value=ypeaks[i];
    m1.insert(make_pair(key,value));//插入key，value值
    mul1.insert(make_pair(key,value));
    mul1.insert(make_pair(key,value));//重复插入key，value值      
}
//查看map中有多少数据
cout<<m1.size()<<","<<mul1.size();
//查找元素--通过迭代器im1进行访问
//  ***im1->first来访问key，使用im->second访问value
//  ***map.find(key)： 查找一个元素，返回键是key的映射的迭代器
Int_t i = 762;
auto im1=m1.find(i);
if(im1!=m1.end()) 
    cout<<im1->first<<", "<<im1->second<<endl;
else 
    cout<<"key="<<i<<" is not found."<<endl;
//正向遍历
//  ***m1.begin()返回指向map头部的迭代器
//  ***m1.end()返回指向map末尾的迭代器
// Hearder--->tail
int num=0;
cout<<"map elements:"<<endl;
for(auto im=m1.begin();im!=m1.end();im++)
{
    if(num>20) break;
    cout<<" "<<im->first<<", "<<im->second<<endl;
    num++;
}

//反向遍历
//  ***m1.rbegin()返回指向map尾部的逆向迭代器
//  ***m1.rend()返回指向map头部的逆向迭代器
// tail --->    Hearder
num=0;
cout<<"map elements:"<<endl;
for(auto im=m1.rbegin();im!=m1.rend();im++)
{
    if(num>20) break;
    cout<<" "<<im->first<<", "<<im->second<<endl;
    num++;
}
//其他成员函数
//  ***maps.clear()清空
//  ***maps.erase()删除一个元素
//  ***maps.lower_bound() 会返回一个迭代器，它指向键值和参数相等或大于参数的第一个元素。
//  ***maps.upper_bound() 也返回一个迭代器，它指向键值大于函数参数的第一个元素。

2.7 在root文件中创建/读取子文件夹

2.7.1 创建

1. 用TDirectoryFIle创建目录结构；

   • 在文件根目录下创建dir1目录：

     fout->cd();//fout为root文件指针，非tree指针哦！！！
     TDirectoryFile * dir1 = new TDirectoryFile("dir1","dir1");
     

   • 在dir1目录下创建子目录：

     dir1->cd();//dir1为目录指针；
     TDirectoryFile * dir1sub = new TDirectoryFile("dir1sub","dir1sub");
     

2. 创建目录完成后，分别在各个目录下（包括直接在fout文件中）声明要存储的对象并分配内存空间（new）；

   • 可以在目录结构中存储TH，TGraph，TTree等：

     TH1I * h0, * h1, * h2, * h1sub;
     fout->cd(); h0 = new TH1I();
     dir1->cd(); h1 = new TH1I();
     dir2->cd(); h2 = new TH1I();
     dir1sub->cd(); h1sub= new TH1I();
     

   • 填图：

     h0->Fill();
     ...
     

   • 写入root文件：

     fout->Write();
     

2.7.2 读取

1. 打开带有目录的root文件：

   TFile * fin = new TFile("test1.root");
   

2. 查看root文件内容：

   fin->ls();
   

3. 得到目录中的TH图的指针：

   TH1I * h0a = (TH1I*)fin->Get("h0");//得到root内的TH
   TH1I * h1a = (TH1I*)fin->Get("dir1/h1");//得到root内dir1目录下的TH
   TH1I * h2a = (TH1I*)fin->Get("dir2/h2");//得到root内dir2目录下的TH
   TH1I * h1suba = (TH1I*)fin->Get("dir1/dir1sub/h1sub");//得到root内dir1的子目录dir1sub下的TH
   

4. 画图：

   h1suba->Draw();
   c1->Draw();
   

2.8 TCut的使用

2.8.1 ROOT终端调用

• 1 View->To0lbar，点击“剪刀”圈出目标成分；
• 2 鼠标放到圈出区域的边缘（一般为黑线），当指针变为小手时，右击，选择SetName（一般为CUTG）；
• 3 画图时在条件参数下加入上步的名字（CUTG），画图即可；

2.8.2 TCut存成脚本，在终端调用

避免查谱时重复cut

• 1 鼠标圈出cut成分，右键存成cut.c文件；

  //r1wvsr1d_cut_noise_300.c
  {
  //========= Macro generated from object: CUTG/Graph
  //========= by ROOT version6.16/00
  
  TCutG *cutg = new TCutG("CUTG",11);  //TCutG类实例化为cutg，CUTG为该cut的名字，11为该cut有11个点
  cutg->SetVarX("r1.d.xhitnec[0]");           //x轴变量名；
  cutg->SetVarY("r1.w.xhitnec[0]");          //y轴变量名；
  cutg->SetTitle("Graph");
  cutg->SetFillStyle(1000);
  cutg->SetPoint(0,402.788,507.273);
  cutg->SetPoint(1,169.913,512.546);
  cutg->SetPoint(2,126.104,291.076);
  cutg->SetPoint(3,246,161.884);
  cutg->SetPoint(4,670.248,143.429);
  cutg->SetPoint(5,882.372,211.979);
  cutg->SetPoint(6,1039.16,272.62);
  cutg->SetPoint(7,794.756,330.624);
  cutg->SetPoint(8,640.274,375.445);
  cutg->SetPoint(9,495.016,454.542);
  cutg->SetPoint(10,402.788,507.273);
  cutg->Draw("");
  }
  //r1wvsr1d_cut_alpha_300.c
  {
  //========= Macro generated from object: CUTG1/Graph
  //========= by ROOT version6.16/00
  
  TCutG *cutg1 = new TCutG("CUTG1",19);    //TCutG类实例化为cutg1，CUTG1为该cut的名字，19为该cut有19个点
  cutg1->SetVarX("r1.d.xhitnec[0]");               //x轴变量名；
  cutg1->SetVarY("r1.w.xhitnec[0]");              //y轴变量名；
  cutg1->SetTitle("Graph");                           //
  cutg1->SetFillStyle(1000);                          //
  cutg1->SetPoint(0,141.046,841.571);        //点1
  cutg1->SetPoint(1,92.962,893.378);          //点2
  cutg1->SetPoint(2,104.983,766.739);
  cutg1->SetPoint(3,329.374,571.025);
  cutg1->SetPoint(4,473.626,499.071);
  cutg1->SetPoint(5,665.961,401.214);
  cutg1->SetPoint(6,866.31,329.261);
  cutg1->SetPoint(7,1238.96,280.332);
  cutg1->SetPoint(8,1707.78,222.769);
  cutg1->SetPoint(9,2100.46,188.231);
  cutg1->SetPoint(10,2192.62,211.256);
  cutg1->SetPoint(11,1703.77,268.819);
  cutg1->SetPoint(12,1218.92,337.895);
  cutg1->SetPoint(13,894.359,406.971);
  cutg1->SetPoint(14,681.989,478.924);
  cutg1->SetPoint(15,437.563,588.294);
  cutg1->SetPoint(16,257.248,755.227);
  cutg1->SetPoint(17,161.081,812.79);
  cutg1->SetPoint(18,141.046,841.571);
  cutg1->Draw("");
  }
  

• 2 将多个cut文件合并，并做适当修改；

{
//========= Macro generated from object: CUTG/Graph
//========= by ROOT version6.16/00
   //r1wvsr1d_cut_alpha_300.c
   TCutG *cutg = new TCutG("CUTG",11);
   cutg->SetVarX("r1.d.xhitnec[0]");
   cutg->SetVarY("r1.w.xhitnec[0]");
   cutg->SetTitle("Graph");
   cutg->SetFillStyle(1000);
   cutg->SetPoint(0,402.788,507.273);
   cutg->SetPoint(1,169.913,512.546);
   cutg->SetPoint(2,126.104,291.076);
   cutg->SetPoint(3,246,161.884);
   cutg->SetPoint(4,670.248,143.429);
   cutg->SetPoint(5,882.372,211.979);
   cutg->SetPoint(6,1039.16,272.62);
   cutg->SetPoint(7,794.756,330.624);
   cutg->SetPoint(8,640.274,375.445);
   cutg->SetPoint(9,495.016,454.542);
   cutg->SetPoint(10,402.788,507.273);
  // cutg->Draw("");        //注释掉
  //}                                //注释掉 
  //{                                //注释掉
  //========= Macro generated from object: CUTG1/Graph
  //========= by ROOT version6.16/00
   //r1wvsr1d_cut_alpha_300.c 
   //TCutG *cutg1 = new TCutG("CUTG1",19); 
   cutg = new TCutG("CUTG1",19);             //注意此处如何修改
   cutg->SetVarX("r1.d.xhitnec[0]");
   cutg->SetVarY("r1.w.xhitnec[0]");
   cutg->SetTitle("Graph");
   cutg->SetFillStyle(1000);
   cutg->SetPoint(0,141.046,841.571);
   cutg->SetPoint(1,92.962,893.378);
   cutg->SetPoint(2,104.983,766.739);
   cutg->SetPoint(3,329.374,571.025);
   cutg->SetPoint(4,473.626,499.071);
   cutg->SetPoint(5,665.961,401.214);
   cutg->SetPoint(6,866.31,329.261);
   cutg->SetPoint(7,1238.96,280.332);
   cutg->SetPoint(8,1707.78,222.769);
   cutg->SetPoint(9,2100.46,188.231);
   cutg->SetPoint(10,2192.62,211.256);
   cutg->SetPoint(11,1703.77,268.819);
   cutg->SetPoint(12,1218.92,337.895);
   cutg->SetPoint(13,894.359,406.971);
   cutg->SetPoint(14,681.989,478.924);
   cutg->SetPoint(15,437.563,588.294);
   cutg->SetPoint(16,257.248,755.227);
   cutg->SetPoint(17,161.081,812.79);
   cutg->SetPoint(18,141.046,841.571);
   //cutg->Draw("");
  }

• 3 利用.L scriptname.c调用该文件，即可使用其中的CUTG与CUTG1；

2.8.3 TCut存成脚本，在脚本中调用，判断是否在区域内；

• 1 应用“TCut存成脚本，在终端调用”方法得到其脚本；
• 2 使用如下命令判断：

#include "TCutG.h"
#inlcude "scriptname.c"  //或者直接添加在代码中

//得到需要判断的参量
Double_t deltaE = we;
Double_t E = de;
//判断是否在cut内
if(!cutg || cutg->IsInside(E,deltaE))   //!cutg: cutg为空，填图，即cutg未设置；cutg->IsInside(x,y)：判断(x,y)是否在cutg的区域内
{
    pGdata[x][y]->SetPoint(n,E,deltaE);//TGraph * graph->SetPoint(n,x,y):将graph的第n个点设置为(x,y);
}

三、画图操作命令（画图、拟合、坐标修改、颜色更换等等）

3.1 脚本画图

////创建指针为c1的画布，c1可根据需要命名；
    TCanvas *c1 = new TCanvas();
////清空画布c1； 
    c1->Clear();
////画tree中histogram，
////    第一个双引号：内容+座标轴范围与bin，alias为该条件别名；
////    第二个双引号：画图的条件、符合、CUT等；
////    第三个双引号：option，如colz等；           
    tree->Draw("td-tu>>htx(500,-20,50)","","");
////执行画图
////脚本模式中每次更改设置修改画图内容，都需要执行此命令更新执行画图；   
////***后面只说功能，不再写此命令***
    c1->Draw();

3.2 得到histogram或TGraph指针（两种方法注意区分）

为了修改图像或者拟合图像，常常需要将root文件存图，得到图像的指针，方便操作

3.2.1 得到Histogram指针

• 从ipf文件(打开的root文件指针)内得到hctof图指针的方法
  • root环境下可直接用hctof->Draw()；
  • 对于数组型的hitsogram，如hpixelpid[][]，会提示变量未声明，只能用如下方法画图；

TH1D * hh = (TH1D *)ipf->Get("hctof");//得到hctof指针，并给hh
hh->Draw();//显示hh文件内的histogram

• 得到tree->Draw()中Histogram别名的指针，通过指针进行进一步操作

tree->Draw("td-tu>>htx(500,-20,50)","","");
TH1D * htx = (TH1D *)gROOT->FindObject("htx");
htx->Draw();

3.2.2 得到TGraph指针

• 从ipf文件(打开的root文件指针)内得到gctof图指针的方法
  • root环境下可直接用gctof->Draw()；
  • 对于数组型的TGraph，如gpixelpid[][]，会提示变量未声明，只能用如下方法画图；

TGraph * gg = (TGraph *)ipf->Get("gpixellpid[15][13]"); //得到gpixellpid[15][13]指针，并给gg
gg->Draw("AP");//画gg图，AP：只画点，不画线

3.2.3 Histogram转TGraph，以及命名

• 将TH图转换成TGraph，并得到TGraph中事件数

  TGraph *gtarget = new TGraph(tree->GetSelectedRows(),tree->GetV2(), tree->GetV1());
  Double_t event_all = gtarget->GetN();
  

• 给TH或TGraph图谱更改名字

  TH1D * th1 = new TH1D();
  th1->SetNameTitle("name_th1","title_th1");
  TGraph *gr1 = new TGraph();
  gr1->SetNameTitle("name_gr1","title_gr1");
  

注意区分hitstogram与TGraph的区别

• histogram：是类似柱状图的东东，横轴是道址（或其他测量量），纵轴是该量的统计量（该指范围内有几个数）；
• TGraph：是一个点图，图中每个点仅有二维坐标(x,y)，无colz；
• TH1D,TH1F,TH1I：类似histogram图；
• TH2D,TH2F,TH2I：在没有colz的情况下类似TGraph图，但TH2的图的每个点是一个概率分布，而非一个点，特别在colz下比较明显；

3.3 叠图

• 方法一：坐标信息由same前的Draw()定死了；

hty->Draw();
htx->Draw("same");

• 方法二：根据每个图的大小，自动调节坐标范围；

THStack *hs = new THStack("hs","test stacked histograms");
hs->Add(htu);
hs->Add(htd);
hs->Draw("nostack");//无此命令会自动归一到第一个图的坐标；
c1->Draw();

3.4 一个画布画多个图

方法1

TH2D * th1;
TH2D * th2;

tree->Draw("r1.w.xhittc[0][0]-r1.w.xtref:r1.w.xhits>>th1(16,-0.5,15.5,4100,-4000,100)","r1.w.xhit==1","colz");
th1 = (TH2D*)gROOT->FindObject("th1");
tree->Draw("r1.w.xhittc[0][0]-r1.w.xtref:r1.w.xhits>>th2(16,-0.5,15.5,4100,-4000,100)","r1.w.xhit==1&&r1.w.xhittc[0][0]>18500&&r1.w.xtref>22000","colz");
th2 = (TH2D*)gROOT->FindObject("th2");

c1->Clear();
c1->Divide(2,1);

c1->cd(1);
gPad->SetLogz();
th1->Draw("colz");
c1->cd(2);
gPad->SetLogz();
th2->Draw("colz");

c1->cd();
c1->Modified();
c1->Update();
c1->Draw();

方法2

• 以下方法在online上可用，但jupyter无法使用还需要调试！

TCanvas *c = new TCanvas("c","Online",100,100,2300,1150);
c->Draw();
c->cd();
TPad *p1 = new TPad("p1","MADC",0,0.33,0.714,1);
p1->Divide(5,2); p1->Draw();
p1->cd(1); m1->Draw("col");
p1->cd(2); m2->Draw("col");
p1->cd(3); m3->Draw("col");
p1->cd(4); m4->Draw("col");
p1->cd(5); a1->Draw("col");
p1->cd(6); a2->Draw("col");
p1->cd(7); a3->Draw("col");
p1->cd(8); a4->Draw("col");
p1->cd(9); a5->Draw("col");
p1->cd(10); a6->Draw("col");
p1->Modified(); p1->Update();

c->cd();
TPad *p2 = new TPad("p2","GDC",0,0,0.714,0.33);
p2->Divide(2,1); p2->Draw();
p2->GetPad(1)->SetRightMargin(0.02);
p2->GetPad(2)->SetLeftMargin(0.02);
p2->GetPad(2)->SetRightMargin(0.02);
p2->cd(1); g1->Draw();
p2->cd(2); g2->Draw();
p2->Modified(); p2->Update();

3.5 填加图例

auto legend = new TLegend(0.7, 0.7, .9, .9);
legend->SetHeader("Qu with different PID","C"); // option "C" allows to center the header
legend->AddEntry(hquall,"Qu all","l");
legend->AddEntry(hqu0,"Q_u for gamma","l");
legend->AddEntry(hqu1,"Q_u for neutron","l");
legend->AddEntry(hqu2,"Q_u for proton","l");
legend->AddEntry(hqu3,"Q_u for pedal","l");
legend->Draw();

3.6 histogram拟合与结果提取

////tree->Draw()中可在bin前加入别名，后续用此指针(htx)画图、拟合、same等；
    tree->Draw("td-tu>>htx(500,-20,50)","","");
////创建拟合函数，也可用root自带的gaus、pol等；
    TF1 *f1 = new TF1("f1","[0]*TMath::Exp(-0.5*((x-[1])/[2])^2)",39.5,43);
////初始化拟合参量
    f1->SetParameter(0,-350);//大致估计范围
    f1->SetParameter(1,41.5);
    f1->SetParameter(2,0.5);
////拟合histgram,Fit("function","option","",range_low,range_high)
    htx->Fit("f1","R");
    htx->Fit("gaus","","",0, 200); //高斯
    htx->Fit("pol1","","",0, 200); //多项式pol1-9；
////提取拟合结果方法
    TF1 *fgaus[2];
    Double_t ped[2],sigma[2];//u,d
    TString sq[2]={"qu","qd"};
    fgaus[0]=hquall->GetFunction("gaus");//得到拟合函数的指针
    fgaus[1]=hqdall->GetFunction("gaus");
    for(int i=0;i<2;i++) 
    {
        ped[i]=fgaus[i]->GetParameter(1);//得到拟合函数的第二个参数
        sigma[i]=fgaus[i]->GetParameter(2);
////TString的格式化输出。用法与printf一致。
        TString ss;
        ss.Form("ped_%s=%.2f, sigma_%s=%.2f",sq[i].Data(),ped[i],sq[i].Data(),sigma[i]); 
        cout<<ss<<endl;
     }

3.7 脚本存储pdf文件

//1. 画图准备
TCanvas * c1 = new TCanvas("c1","c1",600,400);//画图
c1->SetLogz();//图像设置
gStyle->SetPalette(1);//图像设置
gStyle->SetOptStat(0000000);//图像设置
//2. 存图准备
TString pdffilename;//the fileanme of pdf
pdffilename.Form("./checkPID_%04d_%04d.pdf",startrun,stoprun);
TString printpdf; //for save pdf
TString printpdfopenflag = "[";//flag of open pdf
TString printpdfcloseflag = "]";//flag of close pdf
//3. 开始存图
printpdf = pdffilename + printpdfopenflag;
c1->Print(printpdf.Data());
//... ... ...循环此步骤，存入多个图到一个pdf... ... ...
ipt->Draw("l0.w.xhitnec[0]:l0.d.xhitnec[0]>>(1250,0,5000,2000,0,8000)","l0.w.xhit==1&&l0.d.xhit==1","colz");//画图
c1->Update();//更新
c1->Print(pdffilename);//存入pdf
//... ... ...循环此步骤，存入多个图到一个pdf... ... ...
//4. 存图完毕
printpdf = pdffilename + printpdfcloseflag;
c1->Print(printpdf.Data());

3.8 其他功能

• histogram不显示传递误差：histogram name->Sumw2(0);

dtd->Sumw2(0);

• 更改histogram线颜色、粗细、种类

htx->SetLineColor(kGreen);
htx->SetLineColor(kRed);
htx->SetLineColor(kBlue);
htx->SetLineColor(kGray);
htx->SetLineColor(kBlack);

htx->SetLineWidth(1);

htx->SetLineStyle(1);

• 更改TGraph数据点颜色、大小、种类

graph->SetMarkerColor(8); 
graph->SetMarkerSize(1);
graph->SetMarkerStyle(8);

• 将htx的二维图向X轴投影

TH1D * htx1 = htx->ProjectionX("projx of htx");

• 清空c1画布的内容；

c1->Clear();

• 将坐标轴改为log坐标；

c1->SetLogx();  gPad->SetLogx();
c1->SetLogy();  gPad->SetLogy();
c1->SetLogz();  gPad->SetLogz();

• 将坐标轴改为正常坐标；

c1->SetLogx(0); gPad->SetLogx(0);
c1->SetLogy(0); gPad->SetLogy(0);
c1->SetLogz(0); gPad->SetLogz(0);

• 更换sea color颜色；

gStyle->SetPalette(1);

• 设置参数的别名：tree->SetAlias()

TString squa,sqda;
squa.Form("iqu-%f",ped[0]);
sqda.Form("iqd-%f",ped[1]);
tree->SetAlias("qua",squa.Data());
tree->SetAlias("qda",sqda.Data());
TString stcut="itu>0&&itu<4000&&itd>0&&itd<4000";

scut=scut+" && "+stcut;
tree->Draw("itd-itu:log(qua/qda)",scut.Data(),"colz");
c1->SetLogz();
c1->SetLogy(0);
c1->Draw();

• 在TH（TH1D * tdiff）中提取bin相关信息

tdiff->GetNbinxX();//提取TH中bin的个数
tdiff->GetBinContent([num]);//提取TH中num处bin的计数；
tdiff->GetBinLowEdge([num]);//提取TH中num处bin的低边界；
tdiff->Rebin([num]);//在tdiff的基础上重新并道；

• 得到图中计数

tree->GetEntries("GCUT");//得到给定条件下，上一tree->Draw()命令得到的图中的事件数；

3.9 颜色以及功能项

3.9.1 Colors

• colz补充

gStyle->SetPalette(num)	gStyle->SetPalette(num)	gStyle->SetPalette(num)
kDeepSea=51	kGreyScale=52	kDarkBodyRadiator=53
kBlueYellow= 54	**kRainBow=55**	kInvertedDarkBodyRadiator=56
kBird=57	kCubehelix=58	kGreenRedViolet=59
kBlueRedYellow=60	kOcean=61	kColorPrintableOnGrey=62
kAlpine=63	kAquamarine=64	kArmy=65
kAtlantic=66	kAurora=67	kAvocado=68
kBeach=69	kBlackBody=70	kBlueGreenYellow=71
kBrownCyan=72	kCMYK=73	kCandy=74
kCherry=75	kCoffee=76	kDarkRainBow=77
kDarkTerrain=78	kFall=79	kFruitPunch=80
kFuchsia=81	kGreyYellow=82	kGreenBrownTerrain=83
kGreenPink=84	kIsland=85	kLake=86
kLightTemperature=87	kLightTerrain=88	kMint=89
kNeon=90	kPastel=91	kPearl=92
kPigeon=93	kPlum=94	kRedBlue=95
kRose=96	kRust=97	kSandyTerrain=98
kSienna=99	kSolar=100	kSouthWest=101
kStarryNight=102	kSunset=103	kTemperatureMap=104
kThermometer=105	kValentine=106	kVisibleSpectrum=107
kWaterMelon=108	kCool=109	kCopper=110
kGistEarth=111	kViridis=112	kCividis=113

3.9.2 width and style of Line

3.9.3 width and style of Marker

3.9.4 option of tree draw and graph

3.9.5 option of TGraph fit

3.9.6 option of TH1 fit

四、功能代码

4.1 利用TMath::Sort()对能量信号排序

• TMath::Sort())为ROOT软件对c++的std::sort()封装函数，用于数组的排序；
• 数据分析中能量排序有两种情况：
  • 1. Hit结构数据；
  • 1. 非Hit结构数据；

4.1.1 对于Hit结构数据

• Hit结构中所有的数组均与hit相关，如xhitec[xhit], xhits[xhit]……

代码示例（该代码由北京大学韩家兴同学分享）

//非刻度（归一、源刻度）的hit结构为整型数据
void DSSD::Sorthitec(Int_t n, Int_t *x, Int_t *y)//Sort(hit,hitec,hits)
{
  Int_t *index = new Int_t[n];
  TMath::Sort(n,x,index,1);

  Int_t *tmpx = new Int_t[n];
  Int_t *tmpy = new Int_t[n];
  for(Int_t i = 0;i<n;i++)
  {
    tmpx[i] = x[index[i]];
    tmpy[i] = y[index[i]];
  }
  memcpy(x,tmpx,sizeof(Int_t)*n);
  memcpy(y,tmpy,sizeof(Int_t)*n);

  delete index;
  delete tmpx;
  delete tmpy;
}
//归一后或刻度后的hit结构为Double_t
void DSSD::Sorthitne(Int_t n, Double_t *x, Int_t *y)//Sort(hit,hitec,hits)
{
  Int_t *index = new Int_t[n];
  TMath::Sort(n,x,index,1);

  Double_t *tmpx = new Double_t[n];
  Int_t *tmpy = new Int_t[n];
  for(Int_t i = 0;i<n;i++)
  {
    tmpx[i] = x[index[i]];
    tmpy[i] = y[index[i]];
  }
  memcpy(x,tmpx,sizeof(Double_t)*n);
  memcpy(y,tmpy,sizeof(Int_t)*n);

  delete index;
  delete tmpx;
  delete tmpy;
}

4.1.2 对于非Hit结构数据

• 非Hit结构数据通常为多维数组，如adc[x][y]，y代表对应的条，拥有实际意义；
• 排序前要对y进行临时存储，以便在作hit结构时提取正确的条信息；
• 本示例代码调用了hit结构数据排序的代码

代码示例

void DSSD::Sortec(Int_t n, Int_t *x, Int_t *y)//Sort(NStripMax,ec,id)
{
  //对ec数组排序
  //而ec数组的元素顺序与strip对应
  //为了记住strip，添加一个id数组，用于记录strip
  //整体来看就是手动做了一个hit结构，然后用Sorthitec函数排序
  for(Int_t i=0;i<NStripMax;i++)
  {
    y[i] = i;
  }
  Sorthitec(n,x,y);//调用hit结构排序代码
}

4.2 同一粒子在硅条正背面响应的挑选

• 做此操作前，硅条正背面能量信号已经归一对齐或者刻度完成；
• 硅条正背面能量信号最好由大到小排序

代码示例

#include <TROOT.h>
#include <TMath.h>
int pairxy()
{
  //原始数组声明和初始化 
  Double_t xvalue[5] = {4101.5,1263.2,40.8,37.8,37.5};
  Double_t yvalue[5] = {4673.2,59.5,56.5,43.1,37.1};
  Int_t xstrip[5] = {5,6,0,3,2};
  Int_t ystrip[5] = {4,3,5,4,5};
  Int_t xhit = 0;
  Int_t yhit = 0;

  //原始数组hit值
  for(Int_t i = 0; i < 5; i++)
  {
    if(xvalue[i]>0) xhit++;
    if(yvalue[i]>0) yhit++;
  }
  //新的对齐数组
  Double_t xnew[5];
  Double_t ynew[5];
  Int_t xsnew[5];
  Int_t ysnew[5];
  Int_t hitn;

  //新的对齐数组初始化
  for(Int_t i = 0; i < 5; i++)
  {
    //xnew[i] = -1;
    //ynew[i] = -1;
    xsnew[i] = -1;
    ysnew[i] = -1;
    xnew[i]  = TMath::QuietNaN();
    ynew[i]  = TMath::QuietNaN();
  }
  hitn = 0;

  //寻找相等数据
  Int_t tempj = 0;//用于寻找匹配，重新定义起点
  if(xhit>0&&yhit>0)
  {
    for(Int_t i = 0; i < xhit; i++)
    {  
      for(Int_t j = tempj; j< yhit; j++)
      {
      cout<<"比较前：i= "<<i<<";j="<<j<<";tempj="<<tempj<<endl;
        //if(xvalue[i] == yvalue[j])
        if(abs(xvalue[i] - yvalue[j])<50)
        {
          xnew[hitn] = xvalue[i];
          ynew[hitn] = yvalue[j];
          xsnew[hitn] = xstrip[i];
          ysnew[hitn] = ystrip[j];
          hitn++;
          i++;
          tempj = j+1;//注意j++与++j的应用区别
          if(i == 5) j = 5;//防止i=5时，j<5，此时循环依然进行，会读取空数据（接近0的值）
        }
      cout<<"比较后：i= "<<i<<";j="<<j<<";tempj="<<tempj<<endl;
      cout<<"-----------------------------------"<<endl;
      }
    }
  }
  cout<<"原始数组："<<endl; 
  cout<<"xvalue[5] = {"<<xvalue[0]<<", "<<xvalue[1]<<", "<<xvalue[2]<<", "<<xvalue[3]<<", "<<xvalue[4]<<"};"<<"xhit = "<<xhit<<endl;
  cout<<"yvalue[5] = {"<<yvalue[0]<<", "<<yvalue[1]<<", "<<yvalue[2]<<", "<<yvalue[3]<<", "<<yvalue[4]<<"};"<<"yhit = "<<yhit<<endl;
  cout<<"================================"<<endl;
  cout<<"hitn = "<<hitn<<endl;
  cout<<"xnew[5] = "<<xnew[0]<<", "<<xnew[1]<<", "<<xnew[2]<<", "<<xnew[3]<<", "<<xnew[4]<<", "<<endl;
  cout<<"ynew[5] = "<<ynew[0]<<", "<<ynew[1]<<", "<<ynew[2]<<", "<<ynew[3]<<", "<<ynew[4]<<", "<<endl;
  cout<<"xsnew[5] = "<<xsnew[0]<<", "<<xsnew[1]<<", "<<xsnew[2]<<", "<<xsnew[3]<<", "<<xsnew[4]<<", "<<endl;
  cout<<"ysnew[5] = "<<ysnew[0]<<", "<<ysnew[1]<<", "<<ysnew[2]<<", "<<ysnew[3]<<", "<<ysnew[4]<<", "<<endl;


  return 0;
}

4.3 由硅条正背面响应条(xstrip,ystrip)转换成($\theta, \phi$)

• 做此操作前，硅条正背面能量信号已经归一对齐或者刻度完成；
• 硅条正背面能量信号最好由大到小排序；
• 感谢北京大学陈家豪同学的提醒与帮助；

代码示例

#include "TVector3.h"
#include "TMath.h"
#include "TGraph.h"
int main()
{
  //距离，角度
  //W1
  Double_t w1d[9]   = {70.,70.,70.,70.,70.,70.,70.,70.,70.};//F-Cenral,F-Up,F-Down,F-Left,F-Right,B-Rp,B-Donw,B-Left，B-Right
  Double_t w1RX[9]  = {0.,-53.,53., 0.,  0.,-154.,154.,  0.,   0.};//F-Cenral,F-Up,F-Down,F-Left,F-Right,B-Rp,B-Donw,B-Left，B-Right
  Double_t w1RY[9]  = {0.,  0., 0.,53.,-53.,   0.,  0.,127.,-127.};//F-Cenral,F-Up,F-Down,F-Left,F-Right,B-Rp,B-Donw,B-Left，B-Right
  Double_t w1RZ[9]  = {0.,  0., 0., 0.,  0.,   0.,  0.,  0.,   0.};//F-Cenral,F-Up,F-Down,F-Left,F-Right,B-Rp,B-Donw,B-Left，B-Right
  const Double_t w1xcenter = 7.5;
  const Double_t w1ycenter = 7.5;
  const Double_t w1width = 3.1;

  //结果存储
  Double_t w1ptheta = -1.;
  Double_t w1pphi = -1.;
  Double_t w1D = -1.;
  Double_t w1rx = -1.;
  Double_t w1ry = -1.;
  Double_t w1rz = -1.;
  Int_t w1xs= -1;
  Int_t w1ys = -1;

  auto * w1gra = new TGraph;

  TFile * opf = new TFile("phivstheta.root","RECREATE");
  TTree * opt = new TTree("tree","tree structure");

  opt->Branch("w1ptheta",&w1ptheta,"w1ptheta/D");
  opt->Branch("w1pphi",&w1pphi,"w1pphi/D");
  opt->Branch("w1xs",&w1xs,"w1xs/I");
  opt->Branch("w1ys",&w1ys,"w1ys/I");
  opt->Branch("w1D",&w1D,"w1D/D");
  opt->Branch("w1rx",&w1rx,"w1rx/D");
  opt->Branch("w1ry",&w1ry,"w1ry/D");
  opt->Branch("w1rz",&w1rz,"w1rz/D");

  TVector3 v1(0.,0.,0.); //注意时会TVector3!!!
  Int_t w1entries = 0;

  for(Int_t i = 0; i < 9; i++)  //9组望远镜依次循环
  {
    //W1
    for(Int_t xs = 0; xs < 16; xs++)
    {
      for(Int_t ys = 0; ys < 16; ys++)
      {
         v1.SetXYZ((xs-w1xcenter)*w1width,(ys-w1ycenter)*w1width,w1d[i]);//设置硅条本征的xy位置，并将探测器放置于为w1d处
         v1.RotateX(3.1415926*w1RX[i]/180.);  //绕X轴旋装
         v1.RotateY(3.1415926*w1RY[i]/180.);  //绕Y轴旋转
         v1.RotateZ(3.1415926*w1RZ[i]/180.);  //绕Z轴旋转
         w1xs = xs;
         w1ys = ys;
         w1ptheta = v1.Theta()*180./3.1415926; //得到旋转后坐标对应的theta角（弧度-》度）
         w1pphi   = v1.Phi()*180./3.1415926;   //得到旋转后坐标对应的phi角（弧度-》度）

         w1gra->SetPoint(w1entries,w1ptheta,w1pphi);
         w1entries++;

         w1D = w1d[i];
         w1rx = w1RX[i];
         w1ry = w1RY[i];
         w1rz = w1RZ[i];
         opt->Fill();
      }
    }
  std::cout<<"i = "<<i<<std::endl;
  }

  w1gra->SetMarkerStyle(31);
  w1gra->SetMarkerSize(1);
  w1gra->SetLineWidth(0);

  w1gra->SetNameTitle("w1phivstheta","w1phivstheta");

  w1gra->Write();

  opt->Write();
  opf->Close();
  return 0;
  }

4.4 同一望远镜两层硅条的PID匹配

• 问题由来

  • 之前做hit结构时，对同一硅条正背面能量作了大小排序；
  • 因此，两层硅的事件可能不对应，偏离PID（即，但从能量对应上，容易导致错位）；

• 解决办法

  • tracking：即利用之前得到的每个硅条上的(theta,phi)，逐个事件对比，寻找两层硅条(theta,phi)最小的事件作为最后的pid；
  • 前提：需要由theta，phi信息(即，之前要完成单个硅条正面-背面的能量匹配，确定响应的像素)

• 代码示例

  • 从马凯同学拷贝而来；
  • 根据自己的需要添加了一些所需变量，但主体框架以及代码未动；
• 感谢北京大学陈家豪同学、马凯同学提醒与帮助；
• 感谢北京大学杨彪同学的代码；

void Tele0::MatchWD()
{
    //为了标定是否参与匹配，避免循环时重复匹配
    Bool_t found_w[5] = {false, false, false, false, false};
    Bool_t found_d[5] = {false, false, false, false, false};
    //通过角度（theta，phi）判断，W1，BB7事件是否对应
    //声明临时变量，储存角度差值信息
    Double_t diff_theta = 9999., diff_theta_tmp = 0.;
    Double_t diff_phi = 9999., diff_phi_tmp = 0.;
    //确定与w第i个事件对应的d的哪个事件
    Int_t match = -1;
    //循环遍历，寻找与w第i个事件角度差最小的d的某个事件
    for(Int_t i = 0; i < w.hitn; i++)
    {
        if(found_w[i]) continue;//如果第i个w已经匹配过，跳过，找下一个w；
        diff_theta     = 9999.;
        diff_phi       = 9999.;
        match          = -1;
        diff_theta_tmp = 9999.;
        diff_phi_tmp   = 9999.;
        for(Int_t j = 0; j < d.hitn; j++)
        {
            if(found_d[j]) continue;//如果第j个d已经匹配过，跳过，找下一个j；
            //将w，d两个角度作差，给如tmp值，准备后续事件做比较，寻找最小差值，即为匹配
            diff_theta_tmp = TMath::Abs(w.theta[i] - d.theta[j]);//做差
            diff_phi_tmp   = TMath::Abs(w.phi[i] - d.phi[i]);
            if(diff_theta > diff_theta_tmp && diff_theta > diff_theta_tmp)//若之前差值，大于新的差值
            {
                diff_theta = diff_theta_tmp;//将新差值赋值给旧差值
                diff_phi = diff_phi_tmp;
                match = j;//标记
            }
        }
        if(match == -1 || diff_theta > 2 || diff_phi > 2) continue;//若最终没有匹配，或匹配的最终值大于2，跳过
        //找到最小值时
            found_w[i] = true;
        found_d[match] = true;
        //新的hit数 = 旧hit数（i，w相关，match，d相关）
           //we[n] = w.xhitnem[i];
           //de[n] = d.xhitnem[match];
           // se  = s.nem;
           Int_t wxstemp = -1;
           Int_t dxstemp = -1;
           wxstemp = w.xhitns[i];
           dxstemp = d.xhitns[match];
           //20220224
           //一次函数刻度
           we[n] = w.cali2b + w.cali2k * (w.duiqib[wxstemp] + w.duiqik[wxstemp] * w.xhitnec[i]);
           de[n] = d.cali2b + d.cali2k * (d.duiqib[dxstemp] + d.duiqik[dxstemp] * d.xhitnec[match]);
           wt[n] = w.xhitntc[i][0];
           dt[n] = d.xhitntc[match][0];
      wgmulti[n] = w.xhitngmulti[i];
      dgmulti[n] = d.xhitngmulti[match];
        theta[n] = d.theta[match];
          phi[n] = d.phi[match];
          wpx[n] = w.xhitns[i]; 
          wpy[n] = w.yhitns[i];
          dpx[n] = d.xhitns[match];
          dpy[n] = d.yhitns[match];
          n++;
    }
    //if(s.ec > 0) se    = s.cali2b + s.cali2k * s.ec;
    if(s.ec > 0) se    = s.cali2b + s.cali2k * s.ec + s.cali2h * s.ec * s.ec;
    else se = 0.;
}

4.5 激发能谱+效率曲线(放大)

• h1->Rebin(Int_t n);//并道
• h1->Smooth(Double_t a);//平滑
• h1->GetMaximum();//得到最大的bin的counts数
• h1->Scale(Double_t a); //图谱放大（缩小）

void xiplot(Int_t xihe1tele, Int_t xihe2tele, Int_t xirebin, Int_t xismooth2, Int_t xismooth4, Int_t xismooth6)
{
    TCanvas * c1 = new TCanvas();
    TH1D * h1;
    TH1D * h2;
    TH1D * h3;
    TH1D * h4;
    TH1D * h5;
    TH1D * h6;
    gStyle->SetOptStat(kFALSE);

    c1->Clear();
    c1->Divide(1,3);

    c1->cd(1);
    TString ei_0gs  = TString::Format("24Mgei_0gs[%d][%d]", xihe1tele, xihe2tele);
    TString eff_0gs = TString::Format("24Mgei_0gs_eff[%d][%d]", xihe1tele, xihe2tele);
    ipf1->cd();
    //h1->Clear();
    h1 = (TH1D*)ipf1->Get(ei_0gs);
    h1->Rebin(xirebin);
    h1->SetLineColor(1);
    h1->Draw();
    ipf2->cd();
    //h2->Clear();
    h2 = (TH1D*)ipf2->Get(eff_0gs);
    h2->SetLineColor(2);
    h2->Smooth(xismooth2);
    Double_t scale2 = h1->GetMaximum()/h2->GetMaximum();
    //cout<<h1->GetMaximum()<<endl;
    //cout<<h2->GetMaximum()<<endl;
    h2->Scale(scale2);
    h2->Draw("same");


    c1->cd(2);
    TString ei_1st  = TString::Format("24Mgei_1st[%d][%d]", xihe1tele, xihe2tele);
    TString eff_1st = TString::Format("24Mgei_1st_eff[%d][%d]", xihe1tele, xihe2tele);
    ipf1->cd();
    h3 = (TH1D*)ipf1->Get(ei_1st);
    h3->Rebin(xirebin);
    h3->SetLineColor(1);
    h3->Draw();
    ipf2->cd();
    h4 = (TH1D*)ipf2->Get(eff_1st);
    h4->SetLineColor(2);
    h4->Smooth(xismooth4);
    Double_t scale4 = h3->GetMaximum()/h4->GetMaximum();
    h4->Scale(scale4);
    h4->Draw("same");

    c1->cd(3);
    TString ei_2nd  = TString::Format("24Mgei_2nd[%d][%d]", xihe1tele, xihe2tele);
    TString eff_2nd = TString::Format("24Mgei_2nd_eff[%d][%d]", xihe1tele, xihe2tele);
    ipf1->cd();
    h5 = (TH1D*)ipf1->Get(ei_2nd);
    h5->Rebin(xirebin);
    h5->SetLineColor(1);
    h5->Draw();
    ipf2->cd();
    h6 = (TH1D*)ipf2->Get(eff_2nd);
    h6->SetLineColor(2);
    h6->Smooth(xismooth6);
    Double_t scale6 = h5->GetMaximum()/h6->GetMaximum();
    h6->Scale(scale6);
    h6->Draw("same");

    c1->cd();
    c1->Modified();
    c1->Update();
    c1->Draw();
}

4.6 激发能谱本底拟合(TH1D*)->ShowBackground()

• h1->ShowBackground();
• h1->Add(h2,h3,A,B);h1=A*h2+B*h2

void xieiwithbackground(Int_t xihe1tele, Int_t xihe2tele, Int_t xirebin, Int_t back1, Int_t back2, Int_t back3)
{
    TCanvas * c1 = new TCanvas();
    TH1D * h1;
    TH1D * h2;
    TH1D * h3;
    TH1 * hb1;
    TH1 * hb2;
    TH1 * hb3;
    gStyle->SetOptStat(kFALSE);
    c1->Clear();
    c1->Divide(1,3);

    c1->cd(1);
    TString ei_0gs  = TString::Format("24Mgei_0gs[%d][%d]", xihe1tele, xihe2tele);
    ipf1->cd();
    h1 = (TH1D*)ipf1->Get(ei_0gs);
    h1->Rebin(xirebin);
    hb1= h1->ShowBackground(back1);
    h1->SetLineColor(1);
    h1->Add(h1,hb1,1,-1);
    h1->Draw();

    c1->cd(2);
    TString ei_1st  = TString::Format("24Mgei_1st[%d][%d]", xihe1tele, xihe2tele);
    ipf1->cd();
    h2 = (TH1D*)ipf1->Get(ei_1st);
    h2->Rebin(xirebin);
    hb2= h2->ShowBackground(back2);
    h2->SetLineColor(1);
    h2->Add(h2,hb2,1,-1);
    h2->Draw();

    c1->cd(3);
    TString ei_2nd  = TString::Format("24Mgei_2nd[%d][%d]", xihe1tele, xihe2tele);
    ipf1->cd();
    h3 = (TH1D*)ipf1->Get(ei_2nd);
    h3->Rebin(xirebin);
    hb3= h3->ShowBackground(back3);
    h3->SetLineColor(1);
    h3->Add(h3,hb3,1,-1);
    h3->Draw();

    c1->cd();
    c1->Modified();
    c1->Update();
    c1->Draw();   
}

!jupyter nbconvert roottips_xi.ipynb --to html

[NbConvertApp] Converting notebook roottips_xi.ipynb to html
[NbConvertApp] Writing 534467 bytes to roottips_xi.html