主目录: ./tracking
./tracking
- main.cpp - 主程序
- Makefile -编译方法
./tracking/include/
- tracking.h - MakeClass 生成的头文件
- ana.h - 用户分析头文件
./tracking/src/
- tracking.C -MakeClass生成的源文件
- ana.cpp 用户分析源文件#ifndef ana_h
#define ana_h
//ROOT
#include <TH2.h>
//C++
#include <iostream>
//user
#include <tracking.h> //由MakeClass生成，为了继承该函数
using namespace std;
class ana:public tracking //从tracking.h中tracking的类中继承其成员变量和成员函数
{
public:
Int_t run; //run numbe
Int_t NumX;
Int_t NumY;
Int_t WhichX[5];
Int_t WhichY[5];
//Declaration by user
Double_t xx[5],xz[5],yy[5],yz[5];//0-1A,1-1B,2-2A,3-2B,4-3.
Double_t xxfit[5],yyfit[5];//3 x,y, 0-measured, 1- fitted.
Double_t anode[5];//0-1A,1-1B,2-2A,3-2B,4-3
Double_t dx[5],dy[5];//0-1A,1-1B,2-2A,3-2B,4-3
Double_t tx,ty,tz;//target position
Double_t kx,bx,ky,by;
Double_t c2nx,c2ny;//chi2/ndf for xfit,yfit
//构造函数：将传入的tree传递给tracking基类，完成初始化
ana(Int_t run_number, TTree * tree):tracking(tree),run(run_number){}; //run = run_number
virtual ~ana() {};
virtual void SetBranch(TTree *tree);
virtual void TrackInit();
virtual void SetTrace(TH2D *h,Double_t k,Double_t b,Int_t min,Int_t max);
virtual void Analysis(TTree *tree);//分析函数，作用等价于原Loop函数
//Int_t NumX, NumY;//有效探测器数量:0-5
//Int_t WhichX[NumX], WhichY[NumY];//有效探测器是谁(0-1A,1-1B,2-2A,3-2B,4-3)
};
#endif
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TF1.h>
#include <TFitResult.h>
#include <TGraph.h>
#include "ana.h"
using namespace std;
void ana::SetBranch(TTree *tree)
{
//measured pos
tree->Branch("xx",&xx,"xx[5]/D");//1A,1B,2A,2B,3->0,1,2,3,4,5
tree->Branch("xz",&xz,"xz[5]/D");//1A,1B,2A,2B,3->0,1,2,3,4,5
tree->Branch("yy",&yy,"yy[5]/D");//1A,1B,2A,2B,3->0,1,2,3,4,5
tree->Branch("yz",&yz,"yz[5]/D");//1A,1B,2A,2B,3->0,1,2,3,4,5
//xx,yy get by fitted for calulating the efficiency of PPAC
//efficiency = the constants of detectted/that fitted
tree->Branch("xxfit",&xxfit,"xxfit[5]/D");
tree->Branch("yyfit",&yyfit,"yyfit[5]/D");
//anode,coincident with xx and yy, respectively, for trigger events;
tree->Branch("anode",&anode,"anode[5]/D");
//由触发的PPAC个数，以及具体哪个PPAC
tree->Branch("NumX",&NumX,"NumX/I");
tree->Branch("NumY",&NumY,"NumY/I");
tree->Branch("WhichX",&WhichX,"WhichX[NumX]/I");
tree->Branch("WhichY",&WhichY,"WhichY[NumY]/I");
//difference between measured and calculated -for pos resolution.
//dx = xx - xxfit
tree->Branch("dx",&dx,"dx[5]/D");
tree->Branch("dy",&dy,"dy[5]/D");
//target x-y from fitted line at z=0;
tree->Branch("tx",&tx,"tx/D");
tree->Branch("ty",&ty,"ty/D");
tree->Branch("tz",&tz,"tz/D");
//ch2/ndf for linear fitting.
tree->Branch("c2nx",&c2nx,"c2nx/D");
tree->Branch("c2ny",&c2ny,"c2ny/D");
//两种实验获取的触发
tree->Branch("beamTrig",&beamTrig,"beamTrig/I");
tree->Branch("must2Trig",&must2Trig,"must2Trig/I");
//原始文件中的
tree->Branch("targetX",&targetX,"targetX");
tree->Branch("targetY",&targetY,"targetY");
}
void ana::TrackInit()
{
for(Int_t i=0;i<5;i++)
{
xx[i]=PPACF8[i][0];
yy[i]=PPACF8[i][1];
xz[i]=PPACF8[i][2];
yz[i]=PPACF8[i][3];
anode[i]=PPACF8[i][4];
xxfit[i]= -999;
yyfit[i]= -999;
dx[i] = -999;
dy[i] = -999;
WhichX[i] = -1;
WhichY[i] = -1;
}
tx=-999;
ty=-999;
tz=-999;
kx=1.0;
bx=0.0;
ky=1.0;
by=0.0;
//PPAC-x挑出有效（X and anode）PPAC个数以及哪个PPAC
//0->1A,1->1B,2->2A,3->2B,4->3
NumX = 0;
for (Int_t i=0;i<5;i++)
{
if (abs(xx[i])<100.&&anode[i]>-999)
{
WhichX[NumX] = i;
NumX = NumX + 1;
}
}
//排除NumX == 2时，均在PPAC1或均在PPAC2的情况
if(NumX==2 && (WhichX[0]==1 && WhichX[1]==2)) NumX = -2;
if(NumX==2 && (WhichX[0]==3 && WhichX[1]==4)) NumX = -2;
//PPAC-Y挑出有效（Y and anode）PPAC个数以及哪个PPAC
//0->1A,1->1B,2->2A,3->2B,4->3
NumY = 0;
for (Int_t i=0;i<5;i++)
{
if (abs(yy[i])<100.&&anode[i]>-999)
{
WhichY[NumY] = i;
NumY = NumY + 1;
}
}
//排除NumY == 2时，均在PPAC1或均在PPAC2的情况
if(NumY==2 && (WhichY[0]==1 && WhichY[1]==2)) NumY = -2;
if(NumY==2 && (WhichY[0]==3 && WhichY[1]==4)) NumY = -2;
}
void ana::SetTrace(TH2D *h,Double_t k,Double_t b,Int_t min,Int_t max)
{
if(h==0) return;
if(min>=max) return;
for(int i=min;i<max;i++)
{
h->Fill(i,(Int_t)(i*k+b));
}
}
void ana::Analysis(TTree *tree)
{
TH2D *htf8xz=new TH2D("htf8xz","xz trace by ppac",2200,-2000,200,300,-150,150);
TH2D* htf8yz=new TH2D("htf8yz","yz trace by ppac",2200,-2000,200,300,-150,150);
SetBranch(tree);
if (fChain == 0) return;
Long64_t nentries = fChain->GetEntriesFast();
Long64_t nbytes = 0, nb = 0;
for (Long64_t jentry=0; jentry<nentries;jentry++) {
Long64_t ientry = LoadTree(jentry);
if (ientry < 0) break;
nb = fChain->GetEntry(jentry); nbytes += nb;
TrackInit();
//PPAC-X多于两个PPAC触发
if(NumX > 1)
{
//fit x-z trajectory
TFitResultPtr r;
auto * grx = new TGraph;
//传件拟合图像，与NumX有关
for(Int_t i = 0;i<NumX;i++)
{
grx->SetPoint(i,xz[WhichX[i]],xx[WhichX[i]]);
}
//创建一次项拟合函数
TF1 *fx=new TF1("fx","pol1",-2000,0);
//fit option: Q: Quiet mode (minimum printing);
// S: The result of the fit is returned in the TFitResultPtr .
//拟合grx并将结果传递给r
r=grx->Fit(fx,"SQ");
//由拟合曲线得到5个PPAC的拟合位置
//曲线确定，有xz即有xxfit值，故与NumX无关
for(Int_t i = 0;i<5;i++)
{
xxfit[i]=fx->Eval(xz[i]);
}
//由拟合曲线推出零点，即靶的位置
kx = fx->GetParameter(1);
bx = fx->GetParameter(0);
tx = bx/(1+kx);
tz = -tx;
//画径迹图
SetTrace(htf8xz,kx,bx,-1800,0);
//得到有效事件实验值与拟合值的差，与NumX有关,有触发的PPAC赋予新值，无触发的PPAC为初始值-999
for(int i=0;i<NumX;i++)
{
dx[WhichX[i]]=xx[WhichX[i]]-fx->Eval(xz[WhichX[i]]);
}
c2nx=r->Chi2()/r->Ndf();//对于任何程序的自动拟合，原则上都要输出拟合误差进行评估
delete grx;
delete fx;
}
//PPAC-X多于两个PPAC触发
if(NumY > 1)
{
//fit y-z trajectory
TFitResultPtr r;
auto * gry = new TGraph;
//传件拟合图像，与NumY有关
for(Int_t i = 0;i<NumY;i++)
{
gry->SetPoint(i,yz[WhichY[i]],yy[WhichY[i]]);
}
//创建一次项拟合函数
TF1 *fy=new TF1("fy","pol1",-2000,0);
//fit option: Q: Quiet mode (minimum printing);
// S: The result of the fit is returned in the TFitResultPtr .
//拟合grx并将结果传递给r
r=gry->Fit(fy,"SQ");
//由拟合曲线得到5个PPAC的拟合位置
//曲线确定，有xz即有xxfit值，故与NumX无关
for(Int_t i = 0;i<5;i++)
{
yyfit[i]=fy->Eval(yz[i]);
}
//由拟合曲线推出零点，即靶的位置
ky = fy->GetParameter(1);
by = fy->GetParameter(0);
ty = ky * tz + by;
//画径迹图
SetTrace(htf8yz,ky,by,-1800,0);
//得到有效事件实验值与拟合值的差，与NumY有关,有触发的PPAC赋予新值，无触发的PPAC为初始值-999
for(int i=0;i<NumY;i++)
{
dy[WhichY[i]]=yy[WhichY[i]]-fy->Eval(yz[WhichY[i]]);
}
c2ny=r->Chi2()/r->Ndf();//对于任何程序的自动拟合，原则上都要输出拟合误差进行评估
delete gry;
delete fy;
}
tree->Fill();
if(jentry%10000==0) cout<<"processing "<<jentry<<endl;
}
htf8xz->Write();
htf8yz->Write();
tree->Write();
}
#include <iostream>
#include <sstream>
#include <TFile.h>
#include <TTree.h>
#include <TString.h>
#include "ana.h"
using namespace std;
int main(int argc, char** argv)
{
if(argc !=2) {
cout<<"Usage: ./"<<argv[0]<<" run_number "<<endl;
return -1;
}
int run_number = atoi(argv[1]);
TString InputPath, OutputPath, infile, outfile;
InputPath = "../";
OutputPath = "./";
infile.Form("%sf8ppac%03d.root", InputPath.Data(), run_number);
//outfile.Form("%shomework%03d.root", OutputPath.Data(), run_number);
outfile.Form("%shomework2.1.root", OutputPath.Data());
//input
TFile *ipf = new TFile(infile);
if(!ipf->IsOpen()) {
cout<<"Cannot open input file: "<<infile<<endl;
return -1;
}
TTree *ipt = (TTree*)ipf->Get("tree");
//output
TFile *opf = new TFile(outfile,"RECREATE");
TTree *opt = new TTree("tree","ppac tracking");
//注意与method1中main.cpp文件的不同哦！！！
ana *tk=new ana(run_number,ipt);//自己新编写的分析类函数,此处不再是tracking类;
tk->Analysis(opt);//自己编写的分析函数，等价MakeClass生产.C文件的Loop函数；
//
ipf->Close();
opf->Close();
return 1;
}
!jupyter nbconvert homework2.1_Code --to html