KAHIBARO
Discord Login Register

14.4. Derived Quantities

Calculating new observables

In event based analysis you rarely use only the raw variables stored in a TTree. Instead, you construct new observables that better describe the physics you care about or simplify later selections and plots. These new quantities are calculated event by event, inside the same loop where you read the input branches.

A derived quantity is any value computed from existing event variables. Examples in simple analyses might include the sum or difference of energies from two detectors, the magnitude of a 2D position from its $(x, y)$ coordinates, or the ratio of two measured quantities such as $E / p$. In more advanced particle physics analyses you might compute transverse momentum, invariant mass, missing energy, or angular separations, but those topics are covered in dedicated later chapters.

The basic pattern to compute a new observable in a traditional event loop is:

  1. Read the input branches for the current event.
  2. Use C++ expressions to compute the derived variable.
  3. Use the derived variable immediately, for example to fill a histogram or to apply a selection, or store it for later use.

A minimal event loop that calculates a new observable might look like this:

cpp
TTree *tree = (TTree*)file->Get("Events");
float E1, E2;
tree->SetBranchAddress("E1", &E1);
tree->SetBranchAddress("E2", &E2);
TH1F *hESum = new TH1F("hESum", "Total energy;E_{1}+E_{2};Events", 100, 0, 10);
Long64_t nEntries = tree->GetEntries();
for (Long64_t i = 0; i < nEntries; ++i) {
    tree->GetEntry(i);
    float E_sum = E1 + E2;
    hESum->Fill(E_sum);
}

Here E_sum is a new observable that is not stored in the original TTree but is computed for each event from two existing branches.

You can also create new branches in a new TTree or in a snapshot that store these derived observables. This is useful when the same derived quantity is used in many later analyses, or when you want to decouple heavy calculations from light plotting jobs. A simple pattern is:

cpp
TFile *fin  = TFile::Open("input.root");
TTree *tin  = (TTree*)fin->Get("Events");
float E1, E2;
tin->SetBranchAddress("E1", &E1);
tin->SetBranchAddress("E2", &E2);
TFile *fout = TFile::Open("output.root", "RECREATE");
TTree *tout = new TTree("Events", "Events with derived quantities");
float E_sum;
tout->Branch("E_sum", &E_sum, "E_sum/F");
Long64_t nEntries = tin->GetEntries();
for (Long64_t i = 0; i < nEntries; ++i) {
    tin->GetEntry(i);
    E_sum = E1 + E2;
    tout->Fill();
}
fout->Write();
fout->Close();

Now E_sum is a first class variable stored as a branch in a ROOT file, so any later analysis can read it directly.

In modern ROOT, RDataFrame offers a more concise way to define new observables without writing explicit loops. The Define method creates a new column from existing columns:

cpp
ROOT::RDataFrame df("Events", "input.root");
auto df2 = df.Define("E_sum", "E1 + E2");
auto hESum = df2.Histo1D({"hESum", "Total energy;E_{1}+E_{2};Events", 100, 0, 10}, "E_sum");

Here E_sum behaves like a virtual branch. It is computed when needed, for example when filling hESum, but does not require you to manage memory or loops explicitly.

Whenever you design new observables, pay attention to their physical meaning and units. Clearly name them to avoid confusion. For instance, prefer E_sum or E_total over x1, and keep consistent notation such as pT, eta, phi for kinematic variables. Also be aware of domains where formulas are meaningful, for example avoid divisions by zero and invalid arguments to functions such as sqrt or log.

Always check that your derived quantities are defined for all relevant events. Protect against divisions by zero, negative values inside sqrt, or probabilities outside the interval $[0, 1]$. Invalid inputs often lead to nan values and subtle analysis bugs.

Combining branches

Many derived observables are built by combining two or more branches from the TTree. Conceptually, this is simple: you read several related variables for each event and use them together to form new quantities such as sums, differences, products, ratios, angles, or magnitudes.

At the simplest level, combining branches is just writing C++ expressions that involve multiple variables. For example, if your tree stores Cartesian coordinates x and y for each hit, the radial distance from the origin is:

$$
r = \sqrt{x^2 + y^2}.
$$

In an event loop this becomes:

cpp
float x, y;
tree->SetBranchAddress("x", &x);
tree->SetBranchAddress("y", &y);
TH1F *hR = new TH1F("hR", "Radial distance; r; Events", 100, 0, 100);
for (Long64_t i = 0; i < tree->GetEntries(); ++i) {
    tree->GetEntry(i);
    float r = std::sqrt(x * x + y * y);
    hR->Fill(r);
}

Sometimes the combination involves a simple algebraic operation such as sum or difference. Other times it may involve more complex formulas, for example time of flight from two detector times, or an energy asymmetry defined as:

$$
A = \frac{E_1 - E_2}{E_1 + E_2}.
$$

The C++ code is directly parallel to the formula:

cpp
float E1, E2;
tree->SetBranchAddress("E1", &E1);
tree->SetBranchAddress("E2", &E2);
TH1F *hAsym = new TH1F("hAsym", "Energy asymmetry;A;Events", 100, -1, 1);
for (Long64_t i = 0; i < tree->GetEntries(); ++i) {
    tree->GetEntry(i);
    float denom = E1 + E2;
    if (denom == 0) continue;
    float A = (E1 - E2) / denom;
    hAsym->Fill(A);
}

The same idea extends to more complicated structures such as arrays or std::vector branches, which are discussed in detail in later chapters. For example, you might combine all hits in an event by summing over a vector:

cpp
std::vector<float> *E_hits = nullptr;
tree->SetBranchAddress("E_hits", &E_hits);
TH1F *hETotal = new TH1F("hETotal", "Total hit energy;E_{total};Events", 100, 0, 20);
for (Long64_t i = 0; i < tree->GetEntries(); ++i) {
    tree->GetEntry(i);
    float E_total = 0.0;
    for (float E : *E_hits) {
        E_total += E;
    }
    hETotal->Fill(E_total);
}

In RDataFrame, combining branches becomes a matter of declaring the expression in Define. For the asymmetry above you can write:

cpp
ROOT::RDataFrame df("Events", "input.root");
auto df2 = df.Define("A", "((E1 - E2) / (E1 + E2))");
auto hAsym = df2.Filter("E1 + E2 != 0")
                .Histo1D({"hAsym", "Energy asymmetry;A;Events", 100, -1, 1}, "A");

You can chain several Define calls to build observables step by step, which often makes the code easier to read:

cpp
auto df2 = df.Define("E_sum", "E1 + E2")
             .Define("E_diff", "E1 - E2")
             .Define("A", "E_diff / E_sum");

When combining branches, make sure that the variables you use are consistently defined for each event. For example, if some events have missing measurements or sentinel values, you may want to skip them or repair them before using them in derived formulas. Also, sometimes a combination is meaningful only within a particular physical region. In such cases, pair the calculation with a selection cut that restricts the analysis to the valid domain.

Finally, consider storing frequently used combinations as new branches or in a separate processed TTree. This can greatly simplify later analyses. For example, you might create a reduced tree that contains only high level quantities such as total energy, asymmetries, or event categories, and use that tree for fast plotting and fitting.

Views: 13

Comments

Please login to add a comment.

Don't have an account? Register now!