5.1 衰变实验的数据分析I - 基于timestamp的事件重构¶
学习目的：¶
- 基于timestamp的事件重构方法
- 衰变事件的数据分析方法
Experiment:¶
- 188.2MeV $^{40}Ar+^{175}Lu$ fusion-evaporation @ SHANS, IMP
- 测量剩余核的$\alpha$衰变：能量，半衰期，级联关系
实验装置：¶
Detectors：¶
- MWPC+ DSSD(128x48) + SSD(1x3)
- 靶上产生的熔合蒸发反应余核经过SHANS分离，穿过MWPC, 注入到DSSD(stop detector)表面。
- DSSD后面的SSD探测器(veto detector)测量穿过DSSD的束流中的轻粒子，这部分粒子无法被MWPC探测。
- 注入到DSSD上的剩余核产生衰变，DSSD测量带电粒子衰变产物，Si-BOX测量从DSSD逃逸的带电粒子能量。
- SSD和Si-Box探测器数据没有包含在本课件的数据内.
参考文献：¶
- Nuclear Instruments and Methods in Physics Research B 317 (2013) 315–318
DAQ:¶
- 数字化方法记录实验数据
- CAEN V1724 100MHz, 14bit digitizer
- Trigger : 自触发模式，即当某一通道的信号幅度超过设定阈值时，直接记录该通道的能量和时间戳(timestamp)等信息。
Data:¶
alpha.root
- Branch：
- energy /keV
- timestamp /ns
- side: 0-front, 1-back
- strip
mwpc.root
- Branch
- energy /ch
- timestamp /ns
In [1]:
TFile *fds=new TFile("alpha.root");
TTree *tree=(TTree*)fds->Get("tree");
TCanvas *c1=new TCanvas("c1","c1");
Double_t energy;
ULong64_t timestamp;
Int_t side,strip;
tree->SetBranchAddress("energy",&energy);
tree->SetBranchAddress("timestamp",×tamp);
tree->SetBranchAddress("side",&side);
tree->SetBranchAddress("strip",&strip);
tree->Print();
正/背面的能量分布¶
In [2]:
tree->Draw("energy>>(3000,0,30000)","side==0");
c1->Draw();
In [3]:
tree->Draw("energy>>(3000,0,30000)","side==1");
c1->Draw();
事件结构:¶
- 事件格式与之前的数据不同。探测系统的每一个子探测器的信息(如DSSD的某一条的位置(side,strip)，能量(energy)，时间(timestamp))构成一个事件。
- 在后续处理中需要将事件重新组织成常规事件结构。
- 常规事件结构：由每一个入射粒子引发的所有探测器信息
- DSSD: x:strip,energy,timestamp ; y:strip,energy,timestamp;
- MWPC: energy,timestamp
- 常规事件结构：由每一个入射粒子引发的所有探测器信息
- 另外，timestamp不按事件编号的增加而增加，即获取事件的顺序并非是按照事件的发生顺序进行记录的。
- 每个插件内若干个通道共用一个buffer，当某一buffer被填满时，直接将该buffer内数据传输到后端，并清空该buffer。
- 高计数率探测器的数据先传输到后端。
- 需要将事件按照时间顺序进行排序
In [4]:
tree->Scan("side:strip:timestamp:int(energy)","","colsize=15",10,1);
正面和背面事件分别按照时间增加的顺序进行排序，重新组织事件结构。¶
- 分别定义正面和背面信号的map结构mapf1和mapb。
- 将timestamp作为map的key值，其他参数(side,strip,energy)作为value. 当map中插入(key,value)值时，map内按照key的增加顺序(缺省)进行自动排序。
In [5]:
//get timestamp table for front(side=0) and back(side=1) side
struct dssd
{
Double_t energy;
Int_t strip;
};
multimap<ULong64_t, dssd> mapf1, mapb;//map for front and back side
dssd ds;
Long64_t nentries = tree->GetEntriesFast();
for (Long64_t jentry=0; jentry<nentries;jentry++) {
tree->GetEntry(jentry);
ds.energy = energy;
ds.strip = strip;
if(side==0) mapf1.insert(make_pair(timestamp,ds));//key:timestamp, value: ds
if(side==1) mapb.insert(make_pair(timestamp,ds));
if(nentries%100000==0) cout<<jentry<<endl;
}
cout<<"Total number of frontside/backside events = "<<mapf1.size()<<" / "<<mapb.size()<<endl;
fds->Close();
查看事件的时间排序结果¶
In [6]:
TString sout;
auto ifts=mapf1.begin();//map<ULong64_t, dssd>::iterator ifts=;
auto ibts=mapb.begin();
TH2F *hxyraw=new TH2F("hxyraw","front-back correlation for unmatched data",3000,0,30000,3000,0,30000);
sout.Form("%4s%15s%15s%5s%5s%5s%5s"
, "i","front-ts","back-ts"," front-strip", " back-strip", " front-energy", " back-energy");
cout<<sout<<endl;
for(int i=0; i<TMath::Min(mapf1.size(),mapb.size()); i++) {
if(i<20) {
sout.Form("%2d %15llu %15llu %10d %10d %10.1f %10.1f"
, i, ifts->first, ibts->first, ifts->second.strip, ibts->second.strip, ifts->second.energy, ibts->second.energy);
cout<<sout<<endl;
}
hxyraw->Fill(ifts->second.energy,ibts->second.energy);//front-back energy correlation.
ifts++;
ibts++;
}
正面-背面关联¶
- 关联关系不正确
In [7]:
hxyraw->Draw("colz");
c1->Draw();
寻找正面和背面的符合时间窗¶
当多个探测器或一个探测器的不同信号是由一个事件(如某一入射粒子)引发时，它们之间就会有稳定的时间关系(称为关联事件)。
- 即它们之间的时间差集中在固定的范围(如Gauss分布)。
而对于不关联事件，探测器之间的时间差是任意的，因此时间差的分布是均匀分布。
map的lower_bound(),upper_bound()方法
- map::lower_bound(key):返回map中第一个大于或等于key的迭代器指针
- map::upper_bound(key):返回map中第一个大于key的迭代器指针
In [8]:
ULong64_t twindow=100000;//ns
TH1I *hdt=new TH1I("hdt","fts-bts distribution(ns)",20000,-100000,100000);
//time window: [ts-twindow, ts+twindow]
for(auto ifts=mapf1.begin(); ifts!=mapf1.end();ifts++) {
auto ibts1 = mapb.lower_bound(ifts->first-twindow);//left boundary
auto ibts2 = mapb.upper_bound(ifts->first+twindow);//right boundary
for ( ; ibts1 != ibts2; ++ibts1) {
int dt=ifts->first-ibts1->first;
hdt->Fill(dt);
}
}
hdt->Draw();
c1->SetLogy();
c1->Draw();
In [9]:
hdt->GetXaxis()->SetRangeUser(-1000,500);
hdt->Draw();
c1->Draw();
In [10]:
ULong64_t toffset = -550;
multimap<ULong64_t,dssd> mapf;
for(auto ifts=mapf1.begin(); ifts!=mapf1.end();ifts++)
mapf.insert(pair<ULong64_t,dssd>(ifts->first -toffset, ifts->second));
In [11]:
cout<<mapf.size()<<endl;
在符合时间窗内观察正背面关联¶
- 可以看到正确的关联关系
- 正面和背面都有相邻条能量共享
In [12]:
TH2F *hxyc=new TH2F("hxyc","front-back correlation within coincidence window",3000,0,30000,3000,0,30000);
In [13]:
twindow=100;//ns
for(auto ifts=mapf.begin(); ifts!=mapf.end();ifts++) {
auto ibts1=mapb.lower_bound(ifts->first - twindow);
auto ibts2=mapb.upper_bound(ifts->first + twindow);
for ( ; ibts1 != ibts2; ++ibts1) {
int dt=ifts->first - ibts1->first;
hxyc->Fill(ifts->second.energy,ibts1->second.energy);
}
}
hxyc->Draw("colz");
c1->SetLogy(0);
c1->SetLogz();
c1->Draw();
背面相邻条之间事件关联¶
- 从当前事件开始，先在200ns范围内寻找事件(事件已按timestamp递增顺序排序)
- 观察事件多重性分布
In [14]:
Int_t nhit;
TH1I *hnhit=new TH1I("hnhit"," nhit ",5,0,5);
In [15]:
twindow=200;//ns
ifts=mapb.begin();
while(ifts!=mapb.end() ) {
nhit=0;
ULong64_t t0=ifts->first;
while (ifts!=mapb.end() ) {
if(ifts->first > t0+twindow) break;
nhit++;
ifts++;
}
hnhit->Fill(nhit);
}
hnhit->Draw("colz");
gPad->SetLogy();
c1->Draw();
观察背面相邻条关联¶
- 根据多重性分布，选择两重作为相邻条事件。
- 正面的多重性分布和时间差分布于背面类似。
In [16]:
Int_t xhit,yhit;
Int_t xstrip[2],ystrip[2];
ULong64_t xtimestamp[2],ytimestamp[2];
Float_t xenergy[2],yenergy[2];
TH2F *hx2=new TH2F("hx2"," x0-x1 energy correlation",500,0,10000,500,0,10000);
TH1I *hdtx=new TH1I("hdtx"," tx0-tx1 ",20,-150,50);
In [17]:
ifts=mapb.begin();
while(ifts!=mapb.end() ) {
xhit=0;
ULong64_t t0=ifts->first;
while (ifts!=mapb.end() && xhit<2) {
if(ifts->first > t0+twindow) break;
xtimestamp[xhit]= ifts->first;
xstrip[xhit]= ifts->second.strip;
xenergy[xhit]=ifts->second.energy;
xhit++;
ifts++;
}
if(xhit==2) {
hx2->Fill(xenergy[0],xenergy[1]);
Int_t dt=xtimestamp[0]-xtimestamp[1];
hdtx->Fill(dt);
}
}
hx2->Draw("colz");
c1->SetLogz(0);
c1->SetLogy(0);
c1->Draw();
In [18]:
hdtx->Draw();
c1->SetLogy();
c1->Draw();//同一面，不同条之间的时间晃动在 +100-100 范围内
合并正背面相邻条的两重事件，忽略三重事件¶
In [19]:
multimap<ULong64_t,dssd> mapfn,mapbn;
ULong64_t t0;
int mby[5],mfx[5];
In [20]:
twindow=200;//ns
ibts=mapb.begin();
while(ibts!=mapb.end() ) {
yhit=0;
t0=ibts->first;
while (ibts!=mapb.end() && yhit<2) {
if(ibts->first > t0+twindow) break;
ytimestamp[yhit]= ibts->first;
ystrip[yhit]= ibts->second.strip;
yenergy[yhit]=ibts->second.energy;
yhit++;
ibts++;
}
if(yhit==1) {
ds.strip=ystrip[0];
ds.energy=yenergy[0];
mapbn.insert(pair<ULong64_t,dssd>(ytimestamp[0],ds));
mby[1]++;
}
if(yhit==2 && abs(ystrip[0]-ystrip[1])==1) {
ds.energy=yenergy[0]+yenergy[1];
ds.strip=ystrip[1];
t0=ytimestamp[1];
if(yenergy[0]>yenergy[1]) {
ds.strip=ystrip[0];
t0=ytimestamp[0];
}
mapbn.insert(pair<ULong64_t,dssd>(t0,ds));
mby[2]++;
}
if(yhit==2 && abs(ystrip[0]-ystrip[1])>1) {//skip
ds.strip=ystrip[0];
ds.energy=yenergy[0];
//mapbn.insert(pair<ULong64_t,dssd>(ytimestamp[0],ds));
ds.strip=ystrip[1];
ds.energy=yenergy[1];
//mapbn.insert(pair<ULong64_t,dssd>(ytimestamp[1],ds));
mby[3]++;
}
}
In [21]:
mby
In [22]:
twindow=200;//ns
ifts=mapf.begin();
while(ifts!=mapf.end() ) {
xhit=0;
t0=ifts->first;
while (ifts!=mapf.end() && xhit<2) {
if(ifts->first > t0+twindow) break;
xtimestamp[xhit]= ifts->first;
xstrip[xhit]= ifts->second.strip;
xenergy[xhit]=ifts->second.energy;
xhit++;
ifts++;
}
if(xhit==1) {
ds.strip=xstrip[0];
ds.energy=xenergy[0];
mapfn.insert(pair<ULong64_t,dssd>(xtimestamp[0],ds));
mfx[1]++;
}
if(xhit==2 && abs(xstrip[0]-xstrip[1])==1) {
ds.energy=xenergy[0]+xenergy[1];
ds.strip=xstrip[1];
t0=xtimestamp[1];
if(xenergy[0]>xenergy[1]) {
ds.strip=xstrip[0];
t0=xtimestamp[0];
}
mapfn.insert(pair<ULong64_t,dssd>(t0,ds));
mfx[2]++;
}
if(xhit==2 && abs(xstrip[0]-xstrip[1])>1) {
ds.strip=xstrip[0];
ds.energy=xenergy[0];
//mapfn.insert(pair<ULong64_t,dssd>(xtimestamp[0],ds));
ds.strip=xstrip[1];
ds.energy=xenergy[1];
//mapfn.insert(pair<ULong64_t,dssd>(xtimestamp[1],ds));
mfx[3]++;
}
}
In [23]:
mfx
相邻条修正完成后的正背面关联¶
- hxy,对比hxyc, 相邻条修正后正背面关联图，杂散点点明显减少。
- 仍有一些相邻条关联的迹象
- 可能是由于探测器阈值高，损失了部分相邻条的小幅度事件
- 应用正背面关联条件生成hxyfbc
In [24]:
TH2F *hxy=new TH2F("hxy","front-back correlation for corrected data",3000,0,30000,3000,0,30000);
TH2F *hxyfbc=new TH2F("hxyfbc","front-back correlation for corrected data with f-b correlation",3000,0,30000,3000,0,30000);
TH2I *hxystrip=new TH2I("hxystrip","front-back position",128,0,128,48,0,48);
In [25]:
struct dssdxy
{
int xstrip;
int ystrip;
Float_t energy;
};
dssdxy xy;
multimap<ULong64_t, dssdxy> mapdssd;
twindow=200;//ns
for(ifts=mapfn.begin(); ifts!=mapfn.end();ifts++) {
auto ibts1 = mapbn.lower_bound(ifts->first-twindow);//left boundary
auto ibts2 = mapbn.upper_bound(ifts->first+twindow);//right boundary
for ( ; ibts1 != ibts2; ++ibts1) {
int dt=ifts->first-ibts1->first;
hxy->Fill(ifts->second.energy,ibts1->second.energy);//no f-b correlation cut
if(abs(ifts->second.energy-ibts1->second.energy)<500) {
hxyfbc->Fill(ifts->second.energy,ibts1->second.energy);//with f-b correlation cut
xy.xstrip=ifts->second.strip;
xy.ystrip=ibts1->second.strip;
xy.energy=ibts1->second.energy;
mapdssd.insert(pair<ULong64_t, dssdxy>(ibts1->first, xy));
hxystrip->Fill(xy.ystrip,xy.xstrip);
}
}
}
hxy->Draw("colz");
c1->SetLogy(0);
c1->SetLogz();
c1->Draw();
熔合蒸发反应余核的$\alpha-$衰变产物¶
In [26]:
//hxy->GetXaxis()->SetRangeUser(6000,8000);
//hxy->GetYaxis()->SetRangeUser(6000,8000);
hxyfbc->Draw("colz");
c1->Draw();
In [27]:
TH1F *hfe=(TH1F*) hxyfbc->ProjectionX("hfe");
hfe->Draw();
c1->Draw();
In [28]:
//hxystrip->Draw("colz");
//c1->Draw();
MWPC¶
- 按timestamp排序
- 选择MWPC-DSSD符合时间窗
In [29]:
TFile *fmw=new TFile("mwpc.root");
TTree *tmw=(TTree*)fmw->Get("tree");
tmw->SetBranchAddress("energy",&energy);
tmw->SetBranchAddress("timestamp",×tamp);
tmw->Print();
In [30]:
tmw->Scan("timestamp:energy","","colsize=15",10,1);
In [31]:
multimap<ULong64_t, Float_t> mapmw1;//map for mwpc
nentries = tmw->GetEntriesFast();
for (Long64_t jentry=0; jentry<nentries;jentry++) {
tmw->GetEntry(jentry);
mapmw1.insert(pair<ULong64_t,Float_t>(timestamp,energy));
if(nentries%100000==0) cout<<jentry<<endl;
}
cout<<"Total number of mwpc events = "<<mapmw1.size()<<endl;
fmw->Close();
In [32]:
auto imw=mapmw1.begin();//
TH1F *hmwe=new TH1F("hmwe","energy of mwpc",3500,0,35000);
for(int i=0; i<mapmw1.size(); i++) {
if(i<10) {
sout.Form("%2d %15llu %5.1f" , i, imw->first, imw->second);
cout<<sout<<endl;
}
hmwe->Fill(imw->second);
imw++;
}
hmwe->Draw();
c1->SetLogy();
c1->Draw();//trheshold:1000
MWPC 与DSSD符合时间窗¶
In [33]:
twindow=5000;//ns
TH1I *hmdt=new TH1I("hmdt","mts-fts distribution(ns)",1000,-5000,5000);
for(auto idts=mapdssd.begin(); idts!=mapdssd.end();idts++) {
auto imts1=mapmw1.lower_bound(idts->first-twindow);
auto imts2=mapmw1.upper_bound(idts->first+twindow);
for ( ; imts1 != imts2; ++imts1) {
int dt=imts1->first-idts->first;
hmdt->Fill(dt);
}
}
hmdt->Draw();
c1->SetLogy();
c1->Draw();
In [34]:
hmdt->GetXaxis()->SetRangeUser(400,1000);
hmdt->Draw();
c1->Draw();
In [35]:
toffset=650;
multimap<ULong64_t, Float_t> mapmw;
for(auto imts=mapmw1.begin(); imts!=mapmw1.end();imts++)
mapmw.insert(pair<ULong64_t,Float_t>(imts->first - toffset, imts->second));
将DSSD和MWPC事件存入ROOT文件¶
- 不在符合时间窗内的mwpc能量记为 me=-1，符合时间窗内时间me>0
In [36]:
TH1F *hmwec=new TH1F("hmwec","energy of mwpc coincided with DSSD ",3500,0,35000);
TH1F *hdsec=new TH1F("hdsec","energy of DSSD coincided with MWPC ",3500,0,30500);
TH1F *hdsenc=new TH1F("hdsenc","energy of DSSD no coincided with MWPC ",3500,0,30500);
TFile *ff=new TFile("sort1.root","recreate");
TTree *t1=new TTree("tree","sorted events");
Float_t de,me;
Int_t x,y;
ULong64_t tt;
t1->Branch("timestamp",&tt,"timestamp/l");
t1->Branch("me",&me,"me/F");
t1->Branch("de",&de,"de/F");
t1->Branch("xstrip",&x,"xstrip/I");
t1->Branch("ystrip",&y,"ystrip/I");
In [37]:
twindow=200;//ns
for(auto idts=mapdssd.begin(); idts!=mapdssd.end();idts++) {
tt=idts->first;
x=idts->second.xstrip;
y=idts->second.ystrip;
de=idts->second.energy;
me=-1;
auto imts1=mapmw.lower_bound(idts->first-twindow);
auto imts2=mapmw.upper_bound(idts->first+twindow);
for( ; imts1!=imts2; imts1++) {
me=imts1->second;
hmwec->Fill(imts1->second);
hdsec->Fill(idts->second.energy);
}
if(me<0) hdsenc->Fill(idts->second.energy);
t1->Fill();
}
hmwec->Draw();
c1->SetLogy();
c1->Draw();
t1->Write();
ff->Close();
与DSSD符合的MWPC能谱¶
MWPC能谱中超界部分在DSSD符合能谱中消失，意味着这部分是由很重的余核引起的，没有穿出MPWC，直接阻停在内部。
In [38]:
hmwe->Draw();
hmwec->SetLineColor(kRed);
hmwec->Draw("same");
c1->SetLogy();
c1->Draw();
与MWPC符合的DSSD能谱¶
从下面图的特征看，alpha峰即熔合蒸发余核的衰变产物只在无MWPC符合时出现。即alpha衰变的产生事件与束流注入事件不同步，衰变事件的时间应在粒子注入之后。将mwpc作为veto，可选出无入射粒子干扰的衰变谱。 0-2000范围内的DSSD能谱是由低能衰变产物和入射轻粒子(p,d等)，引起的。入射轻粒子在MWPC上能量沉积太小，不产生信号，但一般都有足够的能量穿透注入探测器。这部分由注入探测器后的其他Si探测器进行测量，从DSSD谱中将其veto掉(这部分没有包含在当前的数据中)。
In [39]:
hfe->Draw();//dssd only
hdsec->SetLineColor(kGreen);
hdsec->Draw("same");//coincide with MWPC
c1->SetLogy(0);
c1->Draw();
In [40]:
c1->Clear();
hdsenc->Draw();
c1->Draw();
In [41]:
!jupyter nbconvert 5.1_decay_analysis_I.ipynb --to html