DSSD正背面关联（II）-归一

方法1局限

    1. 由于时正背面关联查看的是像素内的事件，因此统计有可能不够；
    1. pedo附近或者低能区域事件堆积，即使使用ROBUST也会带跑拟合；

方法2： 更充分利用x-y面的关联关系

    1. 在x面和y面分别寻找统计最高的条$ (i_{m},j_{m})$;
    1. 以 $j_{m}$ 为参考条，刻度y面一系列条：$ i_{m−k1}, ... ,i_{m}, .. ,i_{m+k2}$;
    1. 以刻度后的$ i_{m−k1}, ... ,i_{m}, .. ,i_{m+k2}$为参考条，刻度 $ j_{m−k3}, ... , j_{m}, ... ,j_{m+k4}$;
    1. 按照上述步骤向外推进，直到最外边缘的条的正背面关联有足够的统计做拟合。

具体示例见李老师课间或者本人作业（还未完成）。。。

具体步骤

1. 确定每个探测器正背面（非单条）关联的大致范围并存成cut.C文件

gROOT->Macro("cut1.C");

2. 找到每个探测器x和y面中计数最高的条的编号

example:

THStack *hs1 = new THStack("hs1","DSSD1: Counts in strips");
tree->Draw("d1xs>>h1xs(32,0,32)","d1xhit==1 && d1yhit==1");
tree->Draw("d1ys>>h1ys(32,0,32)","d1xhit==1 && d1yhit==1");
h1ys->SetLineColor(2);
hs1->Add(h1xs);
hs1->Add(h1ys);
xs=h1xs->GetMaximumBin();
ys=h1ys->GetMaximumBin();
cout<<xs<<" "<<ys<<endl;
hs1->Draw("nostack");
c1->Draw();

3. 将满足所需条件的x-y关联事件以TGraph形式存入新ROOT文件（如：gdet1.root）

条件如

    1. 在cut.C内；
    1. 相邻条无数据或道址底；
    1. 正背面关联事件在MADC线性好的区域如(200, 8000);

example:

bool bveto="d1x[11]<100 && d1x[13]<100 && d1y[12]<100 && d1y[14]<100";
bool b1213="d1x[12]>200 && d1x[12]<8000 && d1y[13]>200&& d1y[13]<8000";
bool bcut=bveto && b1213;

//选择x-y方向条的范围，以每条计数大于1500为界限
int xrange[2]={13,21};
int yrange[2]={8,16};

TGraph *grx[32][32];
TFile *fg=new TFile("gdet1.root");

for(int ix=0;ix<32;ix++) {
    for(int iy=0;iy<32;iy++) {
        grx[ix][iy]=new TGraph;
        grx[ix][iy]=(TGraph*)fg->Get(Form("g%02d%02d",ix,iy));   
    }
}

grx[10][10]->Draw("A*");
c1->Draw();

4. 以y方向ys条为标准刻度x方向xrange[0]-xrange[1]

example:

//刻度准备
TGraph *grxc[32];
TF1 *fpa1= new TF1("fp1","pol1",200,8000);
TF1 *fpa2= new TF1("fp2","pol2",200,8000);
TH2F *hres[32];
TString shead,sout,hresname;
double par[2],par2[32][3], chi2ndf[32];
int xmax[32],Ncx[32];

//开始刻度
shead.Form("%4s%9s%11s%20s%16s%8s%8s","ix","p0","p1","p2","chi2/ndf","xmax","Nc");
cout<<shead<<endl;
for(int ix=xrange[0]; ix<=xrange[1]; ix++) {
    grxc[ix]=new TGraph;
    int ng=0;
    //第一次拟合得到初步参数缩小拟合范围
    //fit with pol1 first, make new TGraph for the next fitting..
    if(grx[ix][ys]->GetN()>20) {//only fit the graph with enough statistics
        grx[ix][ys]->Fit(fpa1,"Q ROB");
        fpa1->GetParameters(&par[0]);
        for(int k=0; k<grx[ix][ys]->GetN(); k++) {
            double x=grx[ix][ys]->GetX()[k];
            double y=grx[ix][ys]->GetY()[k];
            if (abs(y-(par[0]+par[1]*x))<20) {
                grxc[ix]->SetPoint(ng,x,y);
                ng++;
            }            
        }
    }
    //第二次拟合，精确拟合作为输出刻度参数
    //fit with pol2
    if(grxc[ix]->GetN()>20) {
        grxc[ix]->Fit(fpa2,"Q");
        fpa2->GetParameters(&par2[ix][0]);
        chi2ndf[ix]=fpa2->GetChisquare()/fpa2->GetNDF();
        hresname.Form("hres%02d",ix);
        hres[ix] = new TH2F(hresname,hresname,100,-50,50,1000,200,8000);
        xmax[ix]=0;
        for(int j = 0;j < grx[ix][ys]->GetN(); j++){
            double x=grx[ix][ys]->GetX()[j];
            double y=grx[ix][ys]->GetY()[j];
            if(x>xmax[ix]) xmax[ix]=x;
            x = y-(par2[ix][0]+par2[ix][1]*x+par2[ix][2]*x*x);
            hres[ix]->Fill(x,y);

        }
        Ncx[ix]=grxc[ix]->GetN();
        sout.Form("%4d%9.2f%11.6f%20e%16.2f%8d%8d"  ,ix  ,par2[ix][0]  ,par2[ix][1]  ,par2[ix][2]  ,chi2ndf[ix], xmax[ix],Ncx[ix]);
        cout<<sout<<endl;
    }
}

5. 合并xrange范围内数据，用残差关联图检验刻度效果

example:

//查看残差
TGraph *grsumx=new TGraph;
TH2F *hsumresx = new TH2F("hsumrex","hsumresx",100,-50,50,1000,200,8000);
int ngs=0;
TGraph *grsy[32];//x<->y
hsumresx->Reset();
for(int ix=xrange[0];ix<=xrange[1];ix++) {
    for(int i=0;i<grx[ix][ys]->GetN();i++) {
        double x=grx[ix][ys]->GetX()[i];
        double y=grx[ix][ys]->GetY()[i];
        x = par2[ix][0]+par2[ix][1]*x+par2[ix][2]*x*x;
        grsumx->SetPoint(ngs,x,y);
        hsumresx->Fill(x-y,y);
        ngs++;
    }
}

hsumresx->Draw("colz");
c1->Draw();//seems OK!

//calibrate x strip, merge x strips, swap x and y for fitting
for(int iy=0;iy<32;iy++) {
    int ngy=0;
    grsy[iy]=new TGraph;
    for(int ix=xrange[0];ix<=xrange[1];ix++) {
        for(int i=0;i<grx[ix][iy]->GetN();i++) {
            double x=grx[ix][iy]->GetX()[i];
            double y=grx[ix][iy]->GetY()[i];
            x = par2[ix][0]+par2[ix][1]*x+par2[ix][2]*x*x;
            grsy[iy]->SetPoint(ngy,y,x);//swap x,y for TGraph fit, x=ky+b;
            ngy++;
        }
    }
}

6. 以刻度后的xrange条为参考条，刻度y:0-31

example:

TGraph *grsyc[32];
double par2y[32][3];
TH2F *hy[32], *hresy[32];
TF1 *fpy1[32],*fpy2[32]; 
//TH2F *hsumresx = new TH2F("hsumrex","hsumresx",100,-50,50,1000,200,8000);
//int ngs=0;

shead.Form("%4s%9s%11s%20s%16s%8s%8s","iy","p0","p1","p2","chi2/ndf","ymax","Nc");
cout<<shead<<endl;
for(int iy=0;iy<32;iy++) {
    TString sname;
    sname.Form("p%2dy1",iy);
    fpy1[iy]=new TF1(sname.Data(),"pol1",200,8000);
    fpy1[iy]->SetLineColor(kGray);
    int ng=0;
    grsyc[iy]=new TGraph;
    //fit with pol1 first, make new TGraph for the next fitting..
    if(grsy[iy]->GetN()>20) {
        grsy[iy]->Fit(fpy1[iy],"Q ROB");
        fpy1[iy]->GetParameters(&par[0]);
        for(int k=0; k<grsy[iy]->GetN(); k++) {
            double x=grsy[iy]->GetX()[k];
            double y=grsy[iy]->GetY()[k];
            if (abs(y-(par[0]+par[1]*x))<10) {
                grsyc[iy]->SetPoint(ng,x,y);
                ng++;
            }            
        }
    }
    sname.Form("h%2dy",iy);
    hy[iy]=new TH2F(sname.Data(),sname.Data(),1000,0,8000,1000,0,8000);
    for(int i = 0;i < grsyc[iy]->GetN();++i)  
        hy[iy]->Fill(grsy[iy]->GetX()[i],grsy[iy]->GetY()[i]);
    //fit with pol2
    sname.Form("p%2dy2",iy);
    fpy2[iy]=new TF1(sname.Data(),"pol2",200,8000);
    fpy2[iy]->SetLineColor(kGray);
    if(grsyc[iy]->GetN()>20) {
        grsyc[iy]->Fit(fpy2[iy],"Q");
        fpy2[iy]->GetParameters(&par2y[iy][0]);
        chi2ndf[iy]=fpy2[iy]->GetChisquare()/fpy2[iy]->GetNDF();
        hresname.Form("hresy%02d",iy);
        hresy[iy] = new TH2F(hresname,hresname,100,-50,50,1000,200,8000);
        xmax[iy]=0;
        for(int j = 0;j < grsyc[iy]->GetN(); j++)
            if(grsyc[iy]->GetX()[j]>xmax[iy]) xmax[iy]=grsyc[iy]->GetX()[j];
        for(int j = 0;j < grsy[iy]->GetN(); j++){
            double x=grsy[iy]->GetX()[j];
            double y=grsy[iy]->GetY()[j];
            x = y-(par2y[iy][0]+par2y[iy][1]*x+par2y[iy][2]*x*x);
           hresy[iy]->Fill(x,y);      
        }
        Ncx[iy]=grsyc[iy]->GetN();
        sout.Form("%4d%9.2f%11.6f%20e%16.2f%8d%8d"  ,iy  ,par2y[iy][0]  ,par2y[iy][1]  ,par2y[iy][2]  ,chi2ndf[iy], xmax[iy],Ncx[iy]);
        cout<<sout<<endl;
    }
}

PS:参考条拟合系数p0，p2应接近与零，p1应接近1；

7. 再次回到x条进行刻度，以刻度后的y面条为参考条，刻度x面：0～31

example:

int ymax[32],Ncy[32];
TGraph *grsx[32];//x<->y
TGraph *grsxc[32];
double par2x[32][3];
TH2F *hx[32], *hresx[32];
TF1 *fpx1[32],*fpx2[32];

shead.Form("%4s%9s%11s%20s%16s%8s%8s","ix","p0","p1","p2","chi2/ndf","xmax","Nc");
cout<<shead<<endl;
for(int ix=0;ix<32;ix++) {
    int ngy=0;
    grsx[ix]=new TGraph;
    for(int iy=yrange[0];iy<=yrange[1];iy++) {
        for(int i=0;i<grx[ix][iy]->GetN();i++) {
            double x=grx[ix][iy]->GetX()[i];
            double y=grx[ix][iy]->GetY()[i];
            y = par2y[iy][0]+par2y[iy][1]*y+par2y[iy][2]*y*y;
            grsx[ix]->SetPoint(ngy,x,y);//
            ngy++;
        }
    }
    TString sname;
    sname.Form("p%2dx1",ix);
    fpx1[ix]=new TF1(sname.Data(),"pol1",200,8000);
    fpx1[ix]->SetLineColor(kGray);
    int ng=0;
    grsxc[ix]=new TGraph;
    //fit with pol1 first, make new TGraph for the next fitting..
    if(grsx[ix]->GetN()>20) {
        grsx[ix]->Fit(fpx1[ix],"Q ROB");
        fpx1[ix]->GetParameters(&par[0]);
        for(int k=0; k<grsx[ix]->GetN(); k++) {
            double x=grsx[ix]->GetX()[k];
            double y=grsx[ix]->GetY()[k];
            if (abs(y-(par[0]+par[1]*x))<10) {
                grsxc[ix]->SetPoint(ng,x,y);
                ng++;
            }            
        }
    }
    sname.Form("h%2dx",ix);
    hx[ix]=new TH2F(sname.Data(),sname.Data(),1000,0,8000,1000,0,8000);
    for(int i = 0;i < grsxc[ix]->GetN();++i)  
        hx[ix]->Fill(grsx[ix]->GetX()[i],grsx[ix]->GetY()[i]);
    //fit with pol2
    sname.Form("p%2dx2",ix);
    fpx2[ix]=new TF1(sname.Data(),"pol2",200,8000);
    fpx2[ix]->SetLineColor(kGray);
    if(grsxc[ix]->GetN()>20) {
        grsxc[ix]->Fit(fpx2[ix],"Q");
        fpx2[ix]->GetParameters(&par2x[ix][0]);
        chi2ndf[ix]=fpx2[ix]->GetChisquare()/fpx2[ix]->GetNDF();
        hresname.Form("hresx%02d",ix);
        hresx[ix] = new TH2F(hresname,hresname,100,-50,50,1000,200,8000);
        ymax[ix]=0;
        for(int j = 0;j < grsxc[ix]->GetN(); j++)
            if( grsxc[ix]->GetX()[j]>ymax[ix] ) ymax[ix]=grsxc[ix]->GetX()[j];
        for(int j = 0;j < grsx[ix]->GetN(); j++){
            double x=grsx[ix]->GetX()[j];
            double y=grsx[ix]->GetY()[j];
            x = y-(par2x[ix][0]+par2x[ix][1]*x+par2x[ix][2]*x*x);
           hresx[ix]->Fill(x,y);      
        }
        Ncy[ix]=grsxc[ix]->GetN();
        sout.Form("%4d%9.2f%11.6f%20e%16.2f%8d%8d"  ,ix  ,par2x[ix][0]  ,par2x[ix][1]  ,par2x[ix][2]  ,chi2ndf[ix], ymax[ix],Ncy[ix]);
        cout<<sout<<endl;
    }
}

8. 检查刻度残差

example:

TGraph *grall=new TGraph;
TH2F *hresall = new TH2F("hresall","hresall",100,-50,50,1000,200,8000);
int ngsa=0;

hresall->Reset();
for(int ix=0;ix<32; ix++) {
   for(int iy=0;iy<32;iy++) {
    for(int i=0;i<grx[ix][iy]->GetN();i++) {
        double x=grx[ix][iy]->GetX()[i];
        double y=grx[ix][iy]->GetY()[i];
        x = par2x[ix][0]+par2x[ix][1]*x+par2x[ix][2]*x*x;
        y = par2y[iy][0]+par2y[iy][1]*y+par2y[iy][2]*y*y;
        grall->SetPoint(ngsa,x,y);
        hresall->Fill(x-y,y);
        ngsa++;
    }
   }
}
hresall->Draw("colz");
c1->SetLogz();
c1->Draw();//seems OK!
  • 从残差分布看高能段有些偏。主要的原因是1000以下信号过多，拟合时低能段权重过大。
    • 解决方法，减少参与拟合的低能段事件数目
    • 可先将TGraph转换成TH2，再将TH2转成TGraph(TH2上每个格子上的点，在TGpaph中表达成一个点)，此时每个区域的权重完全相等。
In [1]:
!jupyter nbconvert example3_5.ipynb --to html

[NbConvertApp] Converting notebook example3_5.ipynb to html

[NbConvertApp] Writing 369850 bytes to example3_5.html