3.5 DSSD front-back Correlation(II) - dssd1¶

  • 对x-y的所有条进行能量刻度 (幅度归一化)
  • 学习进行复杂分析的思路和方法

3.4节回顾¶

  • 能量沉积为E的粒子入射到 $(x_i,y_j)$时$A_{x_i}和A_{y_j}$满足线性关系: $A_{x_i} = K_{x_iy_j}A_{y_j} + B_{x_iy_j} $

幅度归一化的目标：将x,y侧所有条相对于参考条$x_i$条进行归一，得到线性系数。¶

方法1¶

  • 固定参考条 $x_i$ 条不变，按照上述方程依次遍历所有的 $y_j$ 条, 得到 y 侧的归一化系数：

  • $A_{x_i} = K_{x_iy_j}A_{x_j} + B_{x_iy_j} \to $ $\lbrace K_{x_i,y_j},B_{x_i,y_j}\rbrace \quad ,y_j=0,1,...N_y $

  • 利用上一步骤得到的系数，对 y 侧所有条进行刻度：$\lbrace A'_{y_j} = K_{x_i,y_j}A_{y_j} + B_{x_i,y_j} \rbrace\quad ,j=0,1,...N_y$。

  • 以刻度后 y 侧某一条 $y_j$ 条作为参考条，对所有 x 侧的条进行刻度 $ A'_{y_j}=K_{y_j,x_i}A_{x_i} + B_{y_j,x_i}$

经上述步骤，DSSD的 $x$ 和 $y$ 的所有条能量都归一化到第$i$参考条

  • 该方法利用的是每个x条和每个y条重叠区域事件(pixel), 前提是入射到每一个$(x_i,y_j)$ pixel的粒子有足够的统计。

方法2¶

  • 更充分利用 x-y 侧的关联关系

  1. 在 x 侧和 y 侧分别寻找统计最高的条$(x_{i_m},y_{j_m})$

  2. 以$x_{i_m}$为参考条，刻度y侧一系列条： $y_{j_{m-k_1}},. ,y_{j_m},.,y_{j_{m+k_2}}$ , 要求这些组合有足够的统计。

  3. 以刻度后的$y_{j_{m-k_1}},. ,y_{j_m},.,y_{j_{m+k_2}}$ 为参考条 (合并成宽度为$k_1+k_2$的条)，刻度 $x_{i_{m-k_3}},.,x_{i_m},.,x_{i_{m+k_4}}$, 要求这些组合有足够的统计。

  4. 按照上述步骤向外推进，直到最外边缘的条的正背面关联有足够的统计做拟合。

下列分析以DSSD1为例，按照方法2进行刻度¶

In [1]:
//%jsroot on
TCanvas *c1 = new TCanvas("c1","c1");
TFile *fin = new TFile("./data/d1xy.root");
TTree *tree =(TTree *)fin->Get("tree");
c1->SetLogz();
In [2]:
tree->Draw("ye:xe>>(8000,0,8000,8000,0,8000)","","colz");
c1->Draw();
In [3]:
tree->Draw("iy:ix>>hxy(32,0,32,32,0,32)","","colz");
c1->Draw();
In [4]:
tree->Draw("xe>>(1000,0,8000)","","colz");
c1->Draw();
  1. 信号在500ch附近过度集中，导致x-y拟合时低能区域点的权重过大。
  2. 为了兼顾动态范围内的不同区域的拟合精度，需要将权重过大的区域点的密度降低。

以下分析没有按照2.进行调整。¶

In [5]:
c1->Clear();
c1->SetWindowSize(900,150);
c1->Divide(3,1);
c1->cd(1); tree->Draw("iy:ix>>hxya1(32,0,32,32,0,32)","xe<300","colz");
c1->cd(2); tree->Draw("iy:ix>>hxya2(32,0,32,32,0,32)","xe>300 && xe<800","colz");
c1->cd(3); tree->Draw("iy:ix>>hxya3(32,0,32,32,0,32)","xe>800","colz");
c1->Draw();

1.1 在x方向 ix=16的条统计最高，因此将ix=16作为参考条。¶

In [6]:
c1->Clear();
c1->SetWindowSize(900, 300);
c1->Divide(2,1);

c1->cd(1);
tree->Draw("xe:ye>>h1(8000,0,8000,8000,0,8000)","ix==16","colz");
c1->cd(2);
tree->Draw("iy>>h2(32,0,32)","ix==16","colz");//iy:8-16 ，These strips have sufficient statistics.
c1->Draw();
  • ix==16的条件下，iy: 8-16有600个以上的统计。因此第一步以xe[16]作为参考条，拟合xe[16]:ye[iy],iy=8-16，得到y方向的刻度系数。
  • 将3.4的拟合步骤整合成函数gFit
    • 输入TGraph，返回拟合参数，chi2/ndf, statistics, residual-plot
In [7]:
void gFit(vector<double> vx, vector<double> vy, double* par, double &chi2ndf, int *sta)
{    
    if(vx.size()<100) return;//to avoid fitting empty TGraph
    double minres = 20;
    TGraph *g1 = new TGraph(vx.size(),&vx[0],&vy[0]); //TGraph Construction with vectors
    TF1 *fpa1 = new TF1("fp1","pol1",200,8000);
    
    double par1[2];
    //robust fitting    
    g1->Fit(fpa1,"Q ROB");
    fpa1->GetParameters(&par1[0]);
    
    //generate new TGraph in which outliers are eliminated.
    vector<double> vx1,vy1;
    for(int i=0; i<vx.size(); i++) {
        if (abs(vy[i]-(par1[0]+par1[1]*vx[i])) < minres) { 
            vx1.push_back(vx[i]);
            vy1.push_back(vy[i]);
        }           
    }
    
    TGraph *g2 = new TGraph(vx1.size(),&vx1[0],&vy1[0]); //TGraph Construction with vectors
    TF1 *fpa2 = new TF1("fp2","pol2",200,8000);    
    //new fit
    g2->Fit(fpa2,"Q");
    fpa2->GetParameters(&par[0]);
    chi2ndf = fpa2->GetChisquare()/fpa2->GetNDF();
    
    sta[0] = *min_element(vy1.begin(),vy1.end());
    sta[1] = *max_element(vy1.begin(),vy1.end());
    sta[2] = vx1.size();
    delete fpa1;
    delete fpa2;
    delete g1;
    delete g2;
   
}

1.2 拟合xe[16] : ye[iy], iy=8-16，得到y方向的刻度系数(刻度ye)¶

In [8]:
int ix,iy;
Int_t xe,ye;
tree->SetBranchAddress("ix",&ix);
tree->SetBranchAddress("iy",&iy);
tree->SetBranchAddress("xe",&xe);
tree->SetBranchAddress("ye",&ye);
In [9]:
vector<double> vye[9],vxe[9];
int ixm = 16;
int iy1 = 8;
int iy2 = 16;
In [10]:
Long64_t nentries = tree->GetEntriesFast();
for (Long64_t jentry = 0; jentry < nentries; jentry++) {
    tree->GetEntry(jentry);
    if(ix == ixm && iy >= iy1 && iy <= iy2) {           
        vye[iy-iy1].push_back(ye);
        vxe[iy-iy1].push_back(xe);
    }

    if (jentry % 10000 == 0 ) cout << ".";
}
cout<<endl;

......................
In [11]:
double parx[32][3], pary[32][3]; //global 
double chi2ndf; 
int sta[3];
TString shead;
In [12]:
shead.Form("%4s%9s%10s%15s%14s%5s%7s%10s","iy","p0","p1","p2","chi2/ndf","ymin","ymax","Ncounts");
cout << shead << endl;
for(int i=0; i<9; i++){
    int iy = i+iy1;
    gFit(vye[i],vxe[i],pary[iy],chi2ndf,sta);//vxe = k* vye + b
    TString pr;
    pr.Form("%4d%9.2f%10.6f%15e%12.2f%7d%7d%8d",
            iy, pary[iy][0], pary[iy][1], pary[iy][2], chi2ndf, sta[0], sta[1], sta[2]);
    cout<<pr<<endl;    
}

  iy       p0        p1             p2      chi2/ndf ymin   ymax   Ncounts
   8     6.39  0.995425  -5.400275e-07       15.48    211   7135     488
   9     9.97  0.965694  -4.026326e-07       18.45    210   7168     512
  10     9.63  0.990557  -4.149655e-07       19.62    210   7532     511
  11    12.77  0.985100  -5.727369e-07       17.52    216   7009     497
  12    11.03  1.002867  -5.214202e-07       16.21    212   7484     563
  13    13.55  0.981965  -3.491393e-07       18.31    214   7048     507
  14     6.83  0.982707  -3.706952e-07       13.86    206   7442     518
  15    16.42  0.991056  -5.016558e-07       16.65    222   6869     492
  16    -4.99  0.982304  -2.302801e-07       21.61    202   7123     497

1.3 刻度ye[iy], iy=8-16，检查x-y关联。¶

In [13]:
TH2F *hxy1 = new TH2F("hxy1","hxy1",8000,0,8000,8000,0,8000);
TH2F *hresy1 = new TH2F("hresy1","hresy1",60,-30,30,800,0,8000);
In [14]:
Long64_t nentries = tree->GetEntriesFast();
for (Long64_t jentry = 0; jentry < nentries; jentry++) {
    tree->GetEntry(jentry);
    if(ix == ixm && iy >= iy1 && iy <= iy2) {           
        double yce = pary[iy][0]+pary[iy][1]*ye+pary[iy][2]*ye*ye;
        hxy1->Fill(yce,xe);
        hresy1->Fill(yce-xe,xe);
    }

    if (jentry % 10000 == 0 ) cout << ".";
}
cout<<endl;

......................
In [15]:
c1->Clear();
c1->SetWindowSize(900, 300);
c1->Divide(2,1);

c1->cd(1);
hxy1->Draw("colz");
c1->cd(2);
hresy1->Draw("colz");
c1->Draw();

2.1 确定下一步拟合的ix范围¶

  • 合并iy:8-16后，ix的所有条的统计都超过600个，因此下一步拟合全部的ix范围
In [16]:
c1->Clear();
c1->Divide(1,1);
tree->Draw("ix>>(32,0,32)","iy>=8 && iy<=16","colz");//iy:8-16 ，These strips have sufficient statistics.
c1->Draw();

2.2. 合并ye[iy], iy=8-16， 进行拟合(刻度xe)¶

In [17]:
vector<double> vxea[32],vyea[32];

Long64_t nentries = tree->GetEntriesFast();
for (Long64_t jentry = 0; jentry < nentries; jentry++) {
    tree->GetEntry(jentry);
    if(iy >= iy1 && iy <= iy2) {           
        double yce = pary[iy][0]+pary[iy][1]*ye+pary[iy][2]*ye*ye;
        vyea[ix].push_back(yce);
        vxea[ix].push_back(xe);
    }

    if (jentry % 10000 == 0 ) cout << ".";
}
cout<<endl;

......................
In [18]:
shead.Form("%4s%9s%10s%15s%14s%5s%7s%10s","ix","p0","p1","p2","chi2/ndf","xmin","xmax","Ncounts");
cout << shead << endl;
for(int i=0; i<32; i++){
    gFit(vxea[i],vyea[i],parx[i],chi2ndf,sta);//vye = k* vxe + b
    TString pr;
    pr.Form("%4d%9.2f%10.6f%15e%12.2f%7d%7d%8d",
            i, parx[i][0], parx[i][1], parx[i][2], chi2ndf, sta[0], sta[1], sta[2]);
    cout<<pr<<endl;    
}

  ix       p0        p1             p2      chi2/ndf xmin   xmax   Ncounts
   0     8.05  1.002932  -1.576983e-07       26.31    209   5062     364
   1     3.48  1.009350  -1.482820e-07       27.39    215   6097     543
   2    -1.48  1.010758  -1.345282e-07       18.76    195   6341     765
   3     3.11  0.994061  -2.826133e-07       17.54    200   5523    1037
   4    -1.05  1.002044   5.947743e-08       21.33    199   5792    1230
   5     3.98  1.020523  -1.546702e-07       20.08    206   5897    1361
   6     0.08  1.015355  -5.303266e-08       23.29    205   7241    1491
   7     3.97  1.012516  -1.561766e-07       22.21    206   6139    1612
   8    -0.69  1.009843  -4.732550e-09       26.22    204   6555    1762
   9     5.59  1.006674   4.464328e-08       24.85    211   6603    1849
  10     1.76  1.022736   3.575231e-08       24.77    204   6192    2009
  11     4.93  1.017132  -2.966935e-10       25.78    205   6614    2075
  12     0.92  0.995566   1.348446e-07       28.43    206   6981    2258
  13     5.86  1.020941   5.603098e-08       27.48    206   6699    2306
  14    -1.31  1.027883   4.288445e-08       25.95    206   7223    2602
  15     7.68  1.002742   5.323300e-08       27.50    208   7218    3101
  16     0.14  0.999849   2.353033e-08       17.66    196   7531    4589
  17     5.78  1.006905  -1.375588e-07       19.23    207   6896    3607
  18    -0.89  1.014668  -6.696057e-08       20.22    194   7159    3855
  19     2.50  0.982619  -8.312678e-08       19.64    197   7125    3646
  20    -1.37  1.011985  -4.759846e-08       18.40    205   6975    3472
  21     2.55  1.010033  -4.654720e-08       20.51    206   6461    2808
  22    -1.93  1.002021  -9.586752e-08       22.06    200   6844    2667
  23     4.93  1.007446  -1.785771e-07       20.81    206   6885    2445
  24     1.78  0.990671  -8.477935e-10       19.71    199   6741    2180
  25     0.52  0.980434  -2.022549e-07       20.59    196   5875    1876
  26    -0.08  0.985974   7.984788e-08       21.16    197   5146    1621
  27     1.80  0.974154  -9.504070e-08       20.89    196   5729    1527
  28     0.61  0.966104   1.205504e-07       20.10    200   5424    1169
  29    -0.41  1.000934  -1.023831e-07       17.89    202   5727     919
  30     0.54  1.006117  -2.390304e-07       20.22    211   5137     549
  31     1.13  0.970377  -7.746834e-08       31.55    210   5698     424
  • ix=16是刻度的参考条，上一步拟合ix的参数应满足 p0$\approx$0, p1$\approx$1.

2.3 刻度xe[ix], ix=0-31，检查x-y关联。¶

In [19]:
TH2F *hxy2 = new TH2F("hxy2","hxy2",8000,0,8000,8000,0,8000);
TH2F *hresy2 = new TH2F("hresy2","hresy2",60,-30,30,800,0,8000);
In [20]:
Long64_t nentries = tree->GetEntriesFast();
for (Long64_t jentry = 0; jentry < nentries; jentry++) {
    tree->GetEntry(jentry);
    if(iy >= iy1 && iy <= iy2) {           
        double yce = pary[iy][0]+pary[iy][1]*ye+pary[iy][2]*ye*ye;
        double xce = parx[ix][0]+parx[ix][1]*xe+parx[ix][2]*xe*xe;        
        hxy2->Fill(yce,xce);
        hresy2->Fill(yce-xce,xce);
    }

    if (jentry % 10000 == 0 ) cout << ".";
}
cout<<endl;

......................
In [21]:
c1->Clear();
c1->SetWindowSize(900, 300);
c1->Divide(2,1);

c1->cd(1);
c1_1->SetLogz();
hxy2->Draw("colz");
c1->cd(2);
c1_2->SetLogz();
hresy2->Draw("colz");
c1->Draw();

3.1 确定下一步拟合的iy范围¶

  • 合并ix:0-31后，iy的所有条的统计都超过600个，因此下一步拟合全部的iy范围
In [22]:
c1->Clear();
c1->Divide(1,1);
tree->Draw("iy>>(32,0,32)","","colz");
c1->Draw();

3.2 合并ye[iy], iy=0-31， 进行拟合(刻度ye)¶

In [23]:
vector<double> vxeb[32],vyeb[32];

Long64_t nentries = tree->GetEntriesFast();
for (Long64_t jentry = 0; jentry < nentries; jentry++) {
    tree->GetEntry(jentry);         
        double xce = parx[ix][0]+parx[ix][1]*xe+parx[ix][2]*xe*xe;
        vyeb[iy].push_back(ye);
        vxeb[iy].push_back(xce);
    if (jentry % 10000 == 0 ) cout << ".";
}
cout<<endl;

......................
In [24]:
shead.Form("%4s%9s%10s%15s%14s%5s%7s%10s","iy","p0","p1","p2","chi2/ndf","ymin","ymax","Ncounts");
cout << shead << endl;
for(int i=0; i<32; i++){
    gFit(vyeb[i],vxeb[i],pary[i],chi2ndf,sta);//vxe = k* vye + b
    TString pr;
    pr.Form("%4d%9.2f%10.6f%15e%12.2f%7d%7d%8d",
            i, pary[i][0], pary[i][1], pary[i][2], chi2ndf, sta[0], sta[1], sta[2]);
    cout<<pr<<endl;    
}

  iy       p0        p1             p2      chi2/ndf ymin   ymax   Ncounts
   0    20.89  0.957691  -3.410491e-07       25.79    210   6340    1768
   1    14.40  0.977566  -4.588626e-07       20.82    207   6679    2578
   2    14.30  0.972767  -3.003339e-07       22.90    209   6559    3526
   3    19.25  0.977489  -3.965019e-07       20.55    206   7418    4423
   4     6.94  0.979184  -2.234616e-07       19.73    202   6852    5232
   5    13.94  0.988205  -3.861101e-07       19.85    205   7297    5678
   6     8.31  0.984323  -3.801819e-07       21.15    206   7349    5996
   7    17.11  0.995165  -5.775258e-07       20.89    209   6936    6473
   8     6.77  0.994894  -4.442756e-07       20.66    200   7135    6917
   9     9.60  0.965925  -4.363193e-07       21.27    204   7168    6937
  10    10.21  0.990275  -3.751683e-07       21.97    206   7532    7073
  11    13.13  0.984726  -5.093436e-07       21.93    212   7237    7020
  12    11.69  1.001960  -3.725290e-07       22.10    207   7484    7365
  13    13.74  0.982161  -3.931464e-07       22.42    200   7139    7191
  14     6.57  0.982686  -3.620009e-07       22.79    200   7442    7298
  15    16.13  0.991095  -4.961962e-07       21.01    206   7218    6907
  16    -4.39  0.981822  -1.497755e-07       22.63    198   7123    7000
  17    10.10  0.970886  -1.759895e-07       20.99    202   7314    6702
  18    -0.56  1.001834  -1.881685e-07       20.07    199   7194    6251
  19     6.34  0.998239  -2.580657e-07       21.94    204   6895    5880
  20    -9.74  0.982470  -1.880519e-09       19.35    196   6780    5358
  21     2.03  0.967963  -4.483838e-08       18.91    196   6852    5341
  22     0.00  0.000000   0.000000e+00       18.91    196   6852    5341
  23     4.72  0.958043  -3.813594e-07       21.86    198   6049    4633
  24     3.73  0.993144  -1.165608e-07       17.63    200   6834    4081
  25    -1.81  0.975686  -1.623837e-07       21.04    198   5809    3832
  26    -1.45  0.994627  -9.622830e-08       19.27    200   6222    3444
  27     6.64  0.983356  -2.815100e-07       20.64    201   5969    2959
  28    -4.89  0.983909   2.625281e-07       21.03    199   5909    2180
  29    -4.19  0.970001   8.237231e-09       20.30    200   6360    1621
  30    -1.23  0.978590  -3.861252e-07       18.15    197   5375    1078
  31    -0.26  0.971015  -3.105138e-07       20.77    201   6225     909

3.3 刻度ye[iy], iy=0-31，检查x-y关联。¶

In [25]:
TH2F *hxy3 = new TH2F("hxy3","hxy3",8000,0,8000,8000,0,8000);
TH2F *hresy3 = new TH2F("hresy3","hresy3",60,-30,30,800,0,8000);
In [26]:
Long64_t nentries = tree->GetEntriesFast();
for (Long64_t jentry = 0; jentry < nentries; jentry++) {
    tree->GetEntry(jentry); 
    
    double yce = pary[iy][0]+pary[iy][1]*ye+pary[iy][2]*ye*ye;
    double xce = parx[ix][0]+parx[ix][1]*xe+parx[ix][2]*xe*xe;        
    hxy3->Fill(yce,xce);
    hresy3->Fill(yce-xce,xce);


    if (jentry % 10000 == 0 ) cout << ".";
}
cout<<endl;

......................
In [27]:
c1->Clear();
c1->SetWindowSize(900, 300);
c1->Divide(2,1);
c1->cd(1);
c1_1->SetLogz();
hxy3->Draw("colz");
c1->cd(2);
c1_2->SetLogz();
hresy3->Draw("colz");
c1->Draw();