DSSD front-back Correlation(II) - dssd1

In the previous tutorial, we introduced the physical basis of front-back energy correlation in a DSSD and defined a detector-internal relative energy scale. We also discussed the most direct calibration approach, namely the pixel-pixel method, in which the relative coefficients are extracted from purified events belonging to individual pixels.

That method is conceptually simple and physically transparent, because the correlation is established directly between the front-side and back-side signals of a specific pixel. However, in practical analysis, the statistical quality of different pixels is often highly uneven. Some pixels contain abundant events and can be fitted reliably, while others suffer from low statistics and cannot provide stable calibration coefficients.

For this reason, a more robust full-detector strategy is needed. In this tutorial, we move from the local pixel-pixel calibration introduced in 3.4 to a more stable strip-normalization method, in which the relative energy scale is propagated strip by strip from the best-populated regions to the full detector.

---

Method 2: Strip Normalization

The basic idea of strip normalization is to avoid relying on one individual pixel at a time. Instead, we begin with the strip that has the highest event statistics, use it to calibrate a statistically strong strip group on the opposite side, and then progressively propagate the relative scale until all strips are covered.

The practical procedure is:

1. identify the strip with the highest event count;
2. use it to normalize a neighboring group on the opposite detector side;
3. combine the normalized group as an effective reference;
4. use that wider reference to normalize the full detector side;
5. return to the opposite side and complete the remaining calibration.

Compared with the pixel-pixel method, this approach is more stable for full-array calibration because each step uses a broader statistical basis.

This is the route used below for DSSD1.

#include <iostream>
#include <fstream>
#include <cmath>

#include "TFile.h"
#include "TTree.h"
#include "TCanvas.h"
#include "TROOT.h"
#include "TStyle.h"
#include "TGraph.h"
#include "TF1.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TRandom.h"
#include "TMath.h"
#include "TString.h"

using namespace std;

// ---------- notebook-global objects ----------
TFile   *fin  = nullptr;
TTree   *tree = nullptr;
TCanvas *c1   = nullptr;

Int_t ix, iy;
Int_t xe, ye;

double parx[32][2], pary[32][2];

// Residual objects kept in memory for direct display
TH2F   *gResiduals[32] = {nullptr};
TCanvas *gResidualCanvas = nullptr;

%%cpp -d
    
void DeleteIfExists(const char* name) {
    TObject *obj = gROOT->FindObject(name);
    if (obj) delete obj;
}

void ClearResiduals() {
    for (int i = 0; i < 32; ++i) {
        if (gResiduals[i]) {
            delete gResiduals[i];
            gResiduals[i] = nullptr;
        }
    }
}

void InitPars() {
    for (int i = 0; i < 32; ++i) {
        parx[i][0] = 0.0; parx[i][1] = 1.0;   // b_x, k_x
        pary[i][0] = 0.0; pary[i][1] = 1.0;   // b_y, k_y
    }
}

inline double Ex(int ix, double xe) {
    return parx[ix][1] * xe + parx[ix][0];
}

inline double Ey(int iy, double ye) {
    return pary[iy][1] * ye + pary[iy][0];
}

bool KeepPoint(double target_raw,
               double center = 500.0,
               double halfWidth = 100.0,
               double keepProb = 0.20)
{
    if (TMath::Abs(target_raw - center) >= halfWidth) return true;
    return gRandom->Rndm() < keepProb;
}

void DrawResiduals(int id1, int id2, const char* cname = "c_residuals")
{
    if (!gResidualCanvas) {
        gResidualCanvas = new TCanvas(cname, cname, 1000, 1400);
    } else {
        gResidualCanvas->SetName(cname);
        gResidualCanvas->SetTitle(cname);
        gResidualCanvas->Clear();
    }

    gResidualCanvas->Divide(4, 8);

    for (int id = id1; id <= id2; ++id) {
        if (!gResiduals[id]) continue;
        gResidualCanvas->cd(id + 1);
        gResiduals[id]->Draw("colz");
    }

    gResidualCanvas->Update();
}

Input Data for Strip Normalization

The input file used in this tutorial is the purified single-pixel event sample prepared in the previous tutorial. These events were selected by applying strict local isolation cuts in order to suppress adjacent-strip charge sharing and preserve only clean front-back correlation data.

The input tree contains four branches,

$$ (ix,\ iy,\ xe,\ ye), $$

where ix and iy are the front-strip and back-strip IDs, and xe, ye are the corresponding raw amplitudes.

These purified events provide the data foundation for the strip-normalization procedure developed below.

DeleteIfExists("c1");
c1 = new TCanvas("c1", "c1", 900, 700);

fin = new TFile("./data/d1xy.root");
tree = (TTree*)fin->Get("tree");

tree->SetBranchAddress("ix", &ix);
tree->SetBranchAddress("iy", &iy);
tree->SetBranchAddress("xe", &xe);
tree->SetBranchAddress("ye", &ye);

InitPars();
gRandom->SetSeed(20260418);

cout << "Entries = " << tree->GetEntries() << endl;

Entries = 210820

Inspect the Global x-y Correlation

Before starting the strip-wise calibration, we first inspect the global front-back correlation of the purified event sample.

A clear diagonal correlation band indicates that the selected events already preserve the underlying front-back energy consistency. This confirms that the purified sample is suitable for constructing a detector-wide relative normalization.

The purpose of this step is not yet to determine calibration coefficients, but to verify the overall quality and usability of the extracted single-pixel data.

DeleteIfExists("hxy_global");
tree->Draw("ye:xe>>hxy_global(8000,0,8000,8000,0,8000)", "", "colz");
c1->SetLogz();
c1->SetGrid(1,1);
c1->Draw();

Reference Strip Selection

In the strip-normalization method, the calibration must start from the most statistically reliable region. Therefore, the first task is to inspect the strip-count distributions and identify the strips with the highest statistics.

For DSSD1, the count distributions show that the front strip with the largest event population is

$$ ix = 16, $$

while the corresponding dense back-strip region lies approximately in

$$ iy = 8 \sim 16. $$

Therefore, the strip ix = 16 is chosen as the initial x-side reference, and the first y-side target group is taken as the central band iy = 8 ... 16.

This choice defines the starting point of the strip-normalization sequence.

DeleteIfExists("hxy_stripmap");
tree->Draw("iy:ix>>hxy_stripmap(32,0,32,32,0,32)", "", "colz");
c1->SetLogz(0);
c1->Draw();

DeleteIfExists("hx");
DeleteIfExists("hy");

tree->Draw("ix>>hx(32,0,32)", "", "goff");
tree->Draw("iy>>hy(32,0,32)", "", "goff");

auto hx = (TH1F*)gROOT->FindObject("hx");
auto hy = (TH1F*)gROOT->FindObject("hy");

DeleteIfExists("c_strip_counts");
auto c_strip_counts = new TCanvas("c_strip_counts", "Strip counts", 900, 700);
hx->SetLineColor(kBlue+1);
hy->SetLineColor(kRed+1);
hx->SetStats(0);
hy->SetStats(0);
hx->Draw("hist");
hy->Draw("hist same");
c_strip_counts->BuildLegend();
c_strip_counts->SetGrid(1,1);
c_strip_counts->Draw();

cout << "Highest x-strip bin = " << hx->GetMaximumBin() - 1 << endl;
cout << "Highest y-strip bin = " << hy->GetMaximumBin() - 1 << endl;

Highest x-strip bin = 16
Highest y-strip bin = 12

DeleteIfExists("hy_ix16");
tree->Draw("iy>>hy_ix16(32,0,32)", "ix==16", "hist");

DeleteIfExists("c_ix16");
auto c_ix16 = new TCanvas("c_ix16", "iy distribution under ix==16", 900, 700);
auto hy_ix16 = (TH1F*)gROOT->FindObject("hy_ix16");
hy_ix16->SetStats(0);
hy_ix16->SetLineColor(kRed+1);
hy_ix16->Draw("hist");
c_ix16->SetGrid(1,1);
c_ix16->Draw();

The comparison above is used to choose the initial reference strip.
For DSSD1, ix = 16 is taken as the initial x-side reference strip, and the first y-side target group is iy = 8 ... 16.

Density Effect in the Dynamic Range

Before performing the fits, it is useful to inspect the amplitude spectra of the purified events.

The spectra show that the event density is strongly concentrated in the range of about 400–700 channels. If all points are used with equal weight, this high-density interval can dominate the regression and slightly bias the fitted coefficients, especially when the normalization is propagated from one strip group to another.

To reduce this statistical overweight, we apply a simple density-control rule to the raw amplitude variable used as the fitting coordinate:

• when fitting y strips, use target_raw = ye;
• when fitting x strips, use target_raw = xe.

This treatment improves the stability of the strip-wise fits over the full dynamic range without changing the physical meaning of the front-back correlation.

DeleteIfExists("hxe");
DeleteIfExists("hye");

tree->Draw("xe>>hxe(100,0,1000)", "", "goff");
tree->Draw("ye>>hye(100,0,1000)", "", "goff");

auto hxe = (TH1F*)gROOT->FindObject("hxe");
auto hye = (TH1F*)gROOT->FindObject("hye");

DeleteIfExists("c_amp");
auto c_amp = new TCanvas("c_amp", "Amplitude spectra", 900, 700);
hxe->SetStats(0);
hye->SetStats(0);
hxe->SetLineColor(kBlue+1);
hye->SetLineColor(kRed+1);
hxe->Draw("hist");
hye->Draw("hist same");
c_amp->SetLogy();
c_amp->SetGrid(1,1);
c_amp->BuildLegend();
c_amp->Draw();

Strip-wise Fit Procedure

To implement the strip-normalization method, we construct a fitting routine that performs the following tasks:

• builds strip-wise correlation graphs from the purified event sample;
• fits each strip correlation with a linear function;
• updates the parameter tables for the x side and y side;
• stores the corresponding residual plots for immediate inspection.

Two parameter tables are maintained throughout the procedure:

• parx[i][0], parx[i][1] for the front strips;
• pary[j][0], pary[j][1] for the back strips.

These parameters define the strip-wise relative transformations

$$ E_{x,i}^{(\mathrm{rel})} = k_{x,i} A_{x,i} + b_{x,i}, \qquad E_{y,j}^{(\mathrm{rel})} = k_{y,j} A_{y,j} + b_{y,j}. $$

The fit routine updates these parameters step by step as the relative scale is propagated across the detector.

void fit(Int_t ix1, Int_t ix2, Int_t iy1, Int_t iy2, TString fitmethod)
{
    Int_t id1 = (fitmethod == "fity") ? iy1 : ix1;
    Int_t id2 = (fitmethod == "fity") ? iy2 : ix2;
    TString sid = (fitmethod == "fity") ? "iy" : "ix";

    TGraph *g[32];
    Int_t   npt[32];

    for (Int_t i = 0; i < 32; ++i) {
        g[i] = new TGraph();
        npt[i] = 0;
    }

    ClearResiduals();

    Long64_t nentries = tree->GetEntries();
    for (Long64_t jentry = 0; jentry < nentries; ++jentry) {
        tree->GetEntry(jentry);

        if (ix < ix1 || ix > ix2 || iy < iy1 || iy > iy2) continue;

        if (fitmethod == "fity") {
            double Eref = Ex(ix, xe);
            double target_raw = ye;
            if (!KeepPoint(target_raw)) continue;
            g[iy]->SetPoint(npt[iy]++, target_raw, Eref);
        } else if (fitmethod == "fitx") {
            double Eref = Ey(iy, ye);
            double target_raw = xe;
            if (!KeepPoint(target_raw)) continue;
            g[ix]->SetPoint(npt[ix]++, target_raw, Eref);
        } else {
            cout << "Unknown fitmethod: " << fitmethod << endl;
            for (Int_t i = 0; i < 32; ++i) delete g[i];
            return;
        }
    }

    cout << Form("%4s%12s%12s%12s%10s",
                 sid.Data(), "b", "k", "chi2/ndf", "Ncounts") << endl;

    for (Int_t id = id1; id <= id2; ++id) {
        if (npt[id] < 20) {
            printf("Warning: too few points for %s-strip %d\n", sid.Data(), id);
            continue;
        }

        g[id]->Fit("pol1", "Q ROB");
        TF1 *f = g[id]->GetFunction("pol1");
        if (!f) continue;

        double b = f->GetParameter(0);
        double k = f->GetParameter(1);
        double c2n = (f->GetNDF() > 0) ? f->GetChisquare() / f->GetNDF() : 0.0;

        if (fitmethod == "fity") {
            pary[id][0] = b;
            pary[id][1] = k;
        } else {
            parx[id][0] = b;
            parx[id][1] = k;
        }

        gResiduals[id] = new TH2F(
            Form("h2res_%s_%02d", sid.Data(), id),
            Form("Residuals %s-strip %d;Residual;Raw amplitude", sid.Data(), id),
            100, -50, 50,
            800, 0, 8000
        );

        double *x = g[id]->GetX();
        double *y = g[id]->GetY();
        for (Int_t i = 0; i < g[id]->GetN(); ++i) {
            double r = y[i] - (b + k * x[i]);
            gResiduals[id]->Fill(r, x[i]);
        }

        cout << Form("%4d%12.3f%12.6f%12.2f%10d",
                     id, b, k, c2n, npt[id]) << endl;
    }

    DrawResiduals(id1, id2, Form("c_res_%s", fitmethod.Data()));

    for (Int_t i = 0; i < 32; ++i) {
        delete g[i];
    }
}

Interpretation of the Two Fit Directions

The strip-normalization procedure alternates between two fitting directions.

fity:Use the already normalized x side as the reference and determine the parameters $(k_{y,j}, b_{y,j})$ for the target y strips.

For each accepted event, the fitted relation is

$$ E_{x,i}^{(\mathrm{rel})} = k_{y,j} A_{y,j} + b_{y,j}. $$

fitx:Use the already normalized y side as the reference and determine the parameters $(k_{x,i}, b_{x,i})$ for the target x strips.

For each accepted event, the fitted relation is

$$ E_{y,j}^{(\mathrm{rel})} = k_{x,i} A_{x,i} + b_{x,i}. $$

By alternating fity and fitx, the relative scale is propagated from the initial high-statistics reference region to the entire detector.

Three-step Normalization Sequence for DSSD1

The following three-step procedure is the practical implementation of Method 2 for DSSD1. Starting from the highest-statistics strip, the relative energy scale is first established in the central region and then extended outward until both detector sides are fully covered.

Step 1: Calibrate the central y-strip band

Use ix = 16 as the initial x-side reference strip and normalize the central back-strip band iy = 8 ... 16.

This step establishes the first reliable normalized region on the y side.

fit(16,16,8,16,"fity");

  iy           b           k    chi2/ndf   Ncounts
   8       9.163    0.992401     3272.70       473
   9      11.815    0.963376     2986.52       496
  10      12.639    0.987991     3102.46       510
  11      15.815    0.981840     4205.15       521
  12      14.121    0.999646     5127.78       587
  13      16.115    0.979757     3740.62       524
  14       8.729    0.980482     2923.32       499
  15      18.706    0.988406     3534.05       494
  16      -2.964    0.980792     2777.09       478

Step 2: Propagate the scale to all x strips

Use the normalized y-strip band iy = 8 ... 16 as the reference and normalize all x strips.

After this step, the full x side is mapped onto a common relative energy scale.

fit(0,31,8,16,"fitx");

  ix           b           k    chi2/ndf   Ncounts
   0      10.364    1.001004     7157.36       498
   1       4.994    1.007995     4894.18       601
   2       0.139    1.009422     4962.89       780
   3       5.417    0.991677     4619.30       977
   4       0.384    1.001152     4051.40      1156
   5       5.822    1.018988     3957.70      1194
   6       1.626    1.013972     3808.66      1361
   7       5.619    1.011041     3949.01      1445
   8       0.897    1.008702     3413.00      1549
   9       6.425    1.006313     3994.41      1739
  10       2.796    1.022540     3759.85      1886
  11       6.414    1.016497     3743.84      1899
  12       2.007    0.995632     4529.19      2186
  13       6.605    1.021009     4373.85      2235
  14      -0.171    1.027828     3924.19      2555
  15       9.094    1.002685     4101.24      3039
  16       0.109    0.999996     3195.20      4574
  17       6.812    1.006059     3472.38      3700
  18      -0.432    1.014300     3382.67      3789
  19       3.190    0.982080     3595.68      3669
  20      -1.060    1.011771     3232.49      3395
  21       3.434    1.009585     3399.96      2623
  22      -0.969    1.001328     3250.49      2476
  23       6.215    1.006309     3192.83      2186
  24       2.722    0.990224     3159.21      2015
  25       2.449    0.978538     3412.90      1645
  26       1.199    0.985310     4074.02      1465
  27       3.583    0.972846     3864.31      1411
  28       0.995    0.966247     3927.16      1114
  29       1.151    0.999710     4578.36       932
  30       2.498    1.004166     4799.13       669
  31       2.913    0.969300     7611.64       623

Step 3: Complete the y-side normalization

Use the fully normalized x side as the reference and normalize all y strips.

After this step, both the x side and the y side are placed on the same detector-internal relative energy scale.

fit(0,31,0,31,"fity");

  iy           b           k    chi2/ndf   Ncounts
   0      23.072    0.955896     5117.92      2145
   1      16.690    0.975289     4565.84      2734
   2      15.992    0.971132     3897.25      3486
   3      21.175    0.975484     4193.28      4142
   4       8.161    0.978042     3674.63      4734
   5      15.865    0.986256     4110.67      5112
   6      10.186    0.982383     3644.45      5566
   7      19.943    0.992019     4225.27      5985
   8       9.217    0.992392     3844.52      6595
   9      12.095    0.963363     3748.67      6875
  10      12.535    0.988047     4165.21      6916
  11      15.934    0.981732     4055.80      6979
  12      13.828    0.999740     4454.71      7242
  13      15.913    0.979871     3687.76      7038
  14       8.813    0.980466     3745.26      7160
  15      18.749    0.988304     4190.12      6541
  16      -2.906    0.980705     3587.92      6524
  17      11.587    0.969743     3644.11      6162
  18       1.030    1.000568     4131.37      5772
  19       8.209    0.996630     3960.71      5235
  20      -9.431    0.982348     3399.13      4900
  21       2.792    0.967595     4706.16      5309
Warning: too few points for iy-strip 22
  23       6.516    0.956050     5049.32      4949
  24       4.870    0.992208     3716.78      3475
  25      -0.325    0.974663     3602.21      3485
  26      -0.239    0.993834     3975.21      3065
  27       8.662    0.981689     3774.22      2642
  28      -4.379    0.984476     4100.75      2161
  29      -3.078    0.969542     4260.31      1746
  30       0.870    0.976613     4513.59      1205
  31       1.880    0.969092     5138.01      1145

Save and Reload the Calibration Parameters

After the strip-normalization procedure is completed, the fitted calibration parameters should be saved for later use.

The saved parameter file contains the relative calibration coefficients for all front strips and back strips. Reloading these parameters makes it possible to apply the same detector-wide normalization to the full data set in subsequent analysis steps.

This also provides the direct input for the event-level front-back energy normalization assignment.

%%cpp -d
    
void saveParameters(const char *fname="cal_dssd1.txt")
{
    ofstream fout(fname);
    for (int i = 0; i < 32; ++i)
        fout << i << " " << parx[i][0] << " " << parx[i][1] << endl;
    for (int i = 0; i < 32; ++i)
        fout << i << " " << pary[i][0] << " " << pary[i][1] << endl;
    fout.close();
}

void readParameters(const char *fname, double px[32][2], double py[32][2])
{
    ifstream fin(fname);
    int id;
    for (int i = 0; i < 32; ++i) fin >> id >> px[i][0] >> px[i][1];
    for (int i = 0; i < 32; ++i) fin >> id >> py[i][0] >> py[i][1];
    fin.close();
}

saveParameters("cal_dssd1.txt");

Global Validation After Normalization

After all strip parameters are determined, the full purified event sample is recalculated using the fitted relative calibration coefficients.

The purpose of this validation is to test whether DSSD1 has been successfully mapped onto a unified relative energy scale. If the strip-normalization procedure is successful, the corrected front-back correlation should be concentrated around the diagonal

$$ E_{x}^{(\mathrm{rel})} \approx E_{y}^{(\mathrm{rel})}, $$ and the residual distribution

$$ E_{y}^{(\mathrm{rel})} - E_{x}^{(\mathrm{rel})} $$

should be centered near zero over the full dynamic range.

This global check provides the final confirmation that the strip-wise propagation has produced a consistent detector-wide normalization.

void calibrateAndPlot(const char* fname="cal_dssd1.txt")
{
    double px[32][2], py[32][2];
    readParameters(fname, px, py);

    DeleteIfExists("hxy_check");
    DeleteIfExists("hdiff_check");

    auto hxy = new TH2F("hxy_check", "cye:cxe;cxe;cye", 400, 0, 8000, 400, 0, 8000);
    auto hdiff = new TH2F("hdiff_check", "cye-cxe;cxe;cye-cxe", 400, 0, 8000, 200, -200, 200);

    Long64_t nentries = tree->GetEntries();
    for (Long64_t jentry = 0; jentry < nentries; ++jentry) {
        tree->GetEntry(jentry);

        double cxe = px[ix][1] * xe + px[ix][0];
        double cye = py[iy][1] * ye + py[iy][0];

        hxy->Fill(cxe, cye);
        hdiff->Fill(cxe, cye - cxe);
    }

    DeleteIfExists("c_check");
    auto c_check = new TCanvas("c_check", "Global check", 800, 400);
    c_check->Divide(2,1);
    c_check->cd(1); hxy->Draw("colz");
    c_check->cd(2); hdiff->Draw("colz");
    c_check->Update();
}

calibrateAndPlot("cal_dssd1.txt");

---

Discussion: Limitations and the Multi-path Approach The sequential method demonstrated above works well because the propagation depth in a DSSD grid is extremely shallow (typically 1–2 steps from the reference), preventing significant error accumulation.

However, its main limitation is the reliance on a manually fixed macroscopic route (e.g., Reference X → Central Y group → All X). In reality, a target strip intersects multiple calibrated strips simultaneously. By forcing a sequential group-by-group calibration, we do not strictly weight each local pixel fit by its statistical quality, thereby under-utilizing the full correlation network of the detector.

For a mathematically rigorous approach that abandons fixed routes, utilizes all available pixel bridges simultaneously, and automatically optimizes the parameters through inverse-variance weighted multi-path fusion, please refer to the supplementary document:

[3.5 Supplementary Note: Multi-path Calibration with Error Propagation]

Assignment

Front-back Energy Normalization on DSSD1, DSSD2, and DSSD3

In this assignment, perform front-back energy normalization for DSSD1, DSSD2, and DSSD3 using the strip-wise relative calibration parameters obtained with the strip-normalization procedure developed in this tutorial.

After the strip parameters are determined, the normalized detector events should be reorganized into a structured hit format for further analysis.

Task

For each detector:

1. apply the front-strip and back-strip normalization parameters to convert raw amplitudes into normalized amplitudes;
2. sort the front-side and back-side signals in descending order of amplitude;
3. apply a threshold of 100 channels;
4. discard all signals below threshold when constructing the hit list;
5. store the normalized event information in a structured format in an output ROOT file.

Required Output Structure

A typical output structure for DSSD1 is

Int_t   x1hit;
Int_t   x1[x1hit];
double  x1e[x1hit];

Int_t   y1hit;
Int_t   y1[y1hit];
double  y1e[y1hit];

double  x1es, y1es;
int     x1m,  y1m;

The same structure should be extended consistently to DSSD2 and DSSD3.

Here:

• x1hit, y1hit are the multiplicities of accepted x-side and y-side hits after threshold selection;
• x1, y1 are the strip IDs of the accepted hits;
• x1e, y1e are the normalized amplitudes of the accepted hits;
• x1es, y1es are the sums of the calibrated amplitudes of all strips on the x side and y side, respectively;
• x1m, y1m are the multiplicities after merging adjacent-strip hits, where adjacent hits are counted as one.