All in One View

Content from Introduction

---

Last updated on 2026-07-10 | Edit this page

Overview

Questions

• How do I access the inclusive kinematics in the simulation output?

Objectives

• Download files for the tutorial.
• Plot basic distributions for x, y, Q2 with different reconstruction methods.

More detailed instructions on streaming and downloading simulation files can be found in the Analysis Tutorial.

Access Simulation from Jefferson Lab xrootd

---

The preferred method for browsing the simulation output is to use xrootd from within the eic-shell. To browse the directory structure and exit, one can run the commands:

BASH

./eic-shell
xrdfs root://dtn-eic.jlab.org
ls /volatile/eic/EPIC/RECO/
exit

Once you’ve located your desired file, you can copy it to your local system using the xrdcp command:

BASH

xrdcp root://dtn-eic.jlab.org//volatile/eic/EPIC/RECO/path-to-file ./

Callout

Note: For simulation campaigns before January 2025, the destination is /work/eic2/EPIC rather than /volatile/eic/EPIC.

Download files for the next step!

---

Let’s start by downloading our files. We will look at two different files from the 26.02.0 reconstruction campaign: a low Q2 and a high Q2 Neutral Current DIS file.

From within the current eic-shell we can grab our files using -

BASH

xrdcp root://dtn-eic.jlab.org//volatile/eic/EPIC/RECO/26.02.0/epic_craterlake/DIS/NC/18x275/minQ2=1/pythia8NCDIS_18x275_minQ2=1_beamEffects_xAngle=-0.025_hiDiv_1.0001.eicrecon.edm4eic.root ./
xrdcp root://dtn-eic.jlab.org//volatile/eic/EPIC/RECO/26.02.0/epic_craterlake/DIS/NC/18x275/minQ2=1000/pythia8NCDIS_18x275_minQ2=1000_beamEffects_xAngle=-0.025_hiDiv_1.0001.eicrecon.edm4eic.root ./

Note that the ./ at the end is the target location to copy to. Change this as desired.

Inspect the InclusiveKinematics branches

---

From here you you can click around the browser to inspect the basic features of the distributions - look for branches titled “InclusiveKinematics*“.

Open the file in ROOT:

BASH

root -l pythia8NCDIS_18x275_minQ2=1_beamEffects_xAngle=-0.025_hiDiv_1.0001.eicrecon.edm4eic.root
TBrowser b

It may be inconvenient to do everything through the TBrowser if you want to compare the distributions, look at multiple files, or to save the histograms. Using your preferred text editor, create a file with the name PlotDistributions.C, and paste the following code:

CPP

void PlotDistributions(TString filename){
  
  std::vector<TString> recon_method = {"Truth", "Electron", "JB", "DA", "Sigma", "ESigma"};

  // Open the file and retrieve the chain
  auto tree = new TChain("events");
  tree->Add(filename);
  
  for (auto method : recon_method){
    
    auto canvas = new TCanvas();
    canvas->Divide(2,2);
    
    TString branch_name;
    canvas->cd(1);
    
    // Draw a histogram for each variable as reconstructed by each method
    branch_name = TString::Format("InclusiveKinematics%s.x",method.Data());
    tree->Draw(branch_name);
    canvas->cd(2);
    branch_name = TString::Format("InclusiveKinematics%s.y",method.Data());
    tree->Draw(branch_name);
    canvas->cd(3);
    branch_name = TString::Format("InclusiveKinematics%s.Q2",method.Data());
    tree->Draw(branch_name);
    canvas->cd(4);
    branch_name = TString::Format("InclusiveKinematics%s.W",method.Data());
    tree->Draw(branch_name);

    branch_name = TString::Format("InclusiveKinematics%s.pdf",method.Data());
    canvas->Print(branch_name); // Write the canvases to a pdf file
  }
}

You can then run this script as

BASH

root -l PlotDistributions.C\(\"pythia8NCDIS_18x275_minQ2=1_beamEffects_xAngle=-0.025_hiDiv_1.0001.eicrecon.edm4eic.root\"\)

replacing the file name with the name of the file that you want to plot.

This script produces histograms of the distributions of the inclusive kinematic variables when reconstructed using several different reconstruction methods. Assuming that there are no major issues in the reconstruction, these distributions should look similar to the true distribution, with some smearing due to resolution effects.

Challenge

Exercise

Run PlotDistributions.C on both the low Q2 and the high Q2 file you downloaded and compare the distributions produced by the different reconstruction methods.

Show me the solution

Each reconstruction method (Electron, JB, DA, Sigma, ESigma) produces its own set of x, y, Q2, and W histograms. They should all broadly resemble the Truth distribution, with method-dependent smearing. Differences between the low and high Q2 files reflect how each method performs in different regions of phase space, which is explored in the following episodes.

Key Points

• Use xrdfs from within the eic-shell to browse available files from simulations campaigns.
• Use xrdcp from within eic-shell to copy files to your local environment.
• Access the reconstructed kinematics using the InclusiveKinematicsXX branches.

Content from Performance Benchmarks

---

Last updated on 2026-07-10 | Edit this page

Overview

Questions

• How do I determine which reconstruction method I should be using?

Objectives

• Produce plots that benchmark the performance of different reconstruction methods.

Using the podio Reader to process simulation files

---

The collections contained in the simulation output often rely on data types made available by edm4hep and edm4eic. These are based on the podio EDM toolkit, which provides its own tools for reading in event data, though approaches using e.g. TTreeReader or RDataFrame are also possible. The data model contains functions that can make key information more accessible. Take the edm4eic:ReconstructedParticle type (see the edm4eic::ReconstructedParticle reference) as an example:

Go into a ROOT prompt (root -l) and create an edm4eic::ReconstructedParticle object

CPP

#include <edm4eic/ReconstructedParticleCollection.h>
edm4eic::ReconstructedParticle rcp

For such an object you can access the tracks or clusters associated with the reconstructed particle as

CPP

rcp.getTracks()
rcp.getClusters()

which would return a list of the associated tracks/clusters. As our rcp was just initialised, the lists are empty - for the objects in the simulation output this won’t be the case.

If you’re not using data frames, you probably do your analysis in an event loop. An event loop with the podio Reader would look somthing like this

CPP

#include "podio/Frame.h"
#include "podio/Reader.h"
#include "edm4eic/ReconstructedParticleCollection.h"

auto reader = podio::makeReader("some_file.root");

for (size_t i = 0; i < reader.getEntries("events"); i++) {
    const auto event = reader.readNextFrame("events");
    auto& reco_collection = event.get<edm4eic::ReconstructedParticleCollection>("ReconstructedParticles");	
    // Your analysis here
}

Below is a full script to produce some resolution benchmark plots using the InclusiveKinematicsXX branches - copy it into a file called BenchmarkReconstruction.C

CPP

// PODIO
#include "podio/Frame.h"
#include "podio/Reader.h"

// DATA MODEL
#include "edm4eic/InclusiveKinematicsCollection.h"

template <class T>
void BinLogX(T *h)
{
   TAxis *axis = h->GetXaxis();
   int bins = axis->GetNbins();

   Axis_t from = TMath::Log10(axis->GetXmin());
   Axis_t to = TMath::Log10(axis->GetXmax());
   Axis_t width = (to - from) / bins;
   Axis_t *new_bins = new Axis_t[bins + 1];

   for (int i = 0; i <= bins; i++) {
     new_bins[i] = TMath::Power(10, from + i * width);
   }
   axis->Set(bins, new_bins);
   delete[] new_bins;
}

void BenchmarkReconstruction(std::string filename, bool bin_log=false) {

  auto reader = podio::makeReader(filename);

  // Declare benchmark histograms
  TH1F *hResoX_electron = new TH1F("hResoX_electron","Electron method;#Deltax/x;Counts",100,-1,1);
  TH1F *hResoX_jb = new TH1F("hResoX_jb","JB method;#Deltax/x;Counts",100,-1,1);
  TH1F *hResoX_da = new TH1F("hResoX_da","Double Angle method;#Deltax/x;Counts",100,-1,1);
  TH1F *hResoX_sigma = new TH1F("hResoX_sigma","#Sigma method;#Deltax/x;Counts",100,-1,1);
  TH1F *hResoX_esigma = new TH1F("hResoX_esigma","e-#Sigma method;#Deltax/x;Counts",100,-1,1);

  TH1F *hResoY_electron = new TH1F("hResoY_electron","Electron method;#Deltay/y;Counts",100,-1,1);
  TH1F *hResoY_jb = new TH1F("hResoY_jb","JB method;#Deltay/y;Counts",100,-1,1);
  TH1F *hResoY_da = new TH1F("hResoY_da","Double Angle method;#Deltay/y;Counts",100,-1,1);
  TH1F *hResoY_sigma = new TH1F("hResoY_sigma","#Sigma method;#Deltay/y;Counts",100,-1,1);
  TH1F *hResoY_esigma = new TH1F("hResoY_esigma","e-#Sigma method;#Deltay/y;Counts",100,-1,1);

  TH1F *hResoQ2_electron = new TH1F("hResoQ2_electron","Electron method;#DeltaQ2/Q2;Counts",100,-1,1);
  TH1F *hResoQ2_jb = new TH1F("hResoQ2_jb","JB method;#DeltaQ2/Q2;Counts",100,-1,1);
  TH1F *hResoQ2_da = new TH1F("hResoQ2_da","Double Angle method;#DeltaQ2/Q2;Counts",100,-1,1);
  TH1F *hResoQ2_sigma = new TH1F("hResoQ2_sigma","#Sigma method;#DeltaQ2/Q2;Counts",100,-1,1);
  TH1F *hResoQ2_esigma = new TH1F("hResoQ2_esigma","e-#Sigma method;#DeltaQ2/Q2;Counts",100,-1,1);

  TH2F *hResoX_2D_electron = new TH2F("hResoX_2D_electron","Electron method;y;#Deltax/x",30,0.001,1,30,-1,1);
  TH2F *hResoX_2D_jb = new TH2F("hResoX_2D_jb","JB method;y;#Deltax/x",30,0.001,1,30,-1,1);
  TH2F *hResoX_2D_da = new TH2F("hResoX_2D_da","Double Angle method;y;#Deltax/x",30,0.001,1,30,-1,1);
  TH2F *hResoX_2D_sigma = new TH2F("hResoX_2D_sigma","#Sigma method;y;#Deltax/x",30,0.001,1,30,-1,1);
  TH2F *hResoX_2D_esigma = new TH2F("hResoX_2D_esigma","e-#Sigma method;y;#Deltax/x",30,0.001,1,30,-1,1);

  TH2F *hResoY_2D_electron = new TH2F("hResoY_2D_electron","Electron method;y;#Deltay/y",30,0.001,1,30,-1,1);
  TH2F *hResoY_2D_jb = new TH2F("hResoY_2D_jb","JB method;y;#Deltay/y",30,0.001,1,30,-1,1);
  TH2F *hResoY_2D_da = new TH2F("hResoY_2D_da","Double Angle method;y;#Deltay/y",30,0.001,1,30,-1,1);
  TH2F *hResoY_2D_sigma = new TH2F("hResoY_2D_sigma","#Sigma method;y;#Deltay/y",30,0.001,1,30,-1,1);
  TH2F *hResoY_2D_esigma = new TH2F("hResoY_2D_esigma","e-#Sigma method;y;#Deltay/y",30,0.001,1,30,-1,1);

  TH2F *hResoQ2_2D_electron = new TH2F("hResoQ2_2D_electron","Electron method;y;#DeltaQ2/Q2",30,0.001,1,30,-1,1);
  TH2F *hResoQ2_2D_jb = new TH2F("hResoQ2_2D_jb","JB method;y;#DeltaQ2/Q2",30,0.001,1,30,-1,1);
  TH2F *hResoQ2_2D_da = new TH2F("hResoQ2_2D_da","Double Angle method;y;#DeltaQ2/Q2",30,0.001,1,30,-1,1);
  TH2F *hResoQ2_2D_sigma = new TH2F("hResoQ2_2D_sigma","#Sigma method;y;#DeltaQ2/Q2",30,0.001,1,30,-1,1);
  TH2F *hResoQ2_2D_esigma = new TH2F("hResoQ2_2D_esigma","e-#Sigma method;y;#DeltaQ2/Q2",30,0.001,1,30,-1,1);

  // Logarithmic binning on x axis of 2D plots
  if (bin_log){
    BinLogX(hResoX_2D_electron);
    BinLogX(hResoX_2D_jb);
    BinLogX(hResoX_2D_da);
    BinLogX(hResoX_2D_sigma);
    BinLogX(hResoX_2D_esigma);
    
    BinLogX(hResoY_2D_electron);
    BinLogX(hResoY_2D_jb);
    BinLogX(hResoY_2D_da);
    BinLogX(hResoY_2D_sigma);
    BinLogX(hResoY_2D_esigma);
    
    BinLogX(hResoQ2_2D_electron);
    BinLogX(hResoQ2_2D_jb);
    BinLogX(hResoQ2_2D_da);
    BinLogX(hResoQ2_2D_sigma);
    BinLogX(hResoQ2_2D_esigma);
  }
  
  Float_t x_truth, x_electron, x_jb, x_da, x_sigma, x_esigma;
  Float_t y_truth, y_electron, y_jb, y_da, y_sigma, y_esigma;
  Float_t Q2_truth, Q2_electron, Q2_jb, Q2_da, Q2_sigma, Q2_esigma;

  cout << reader.getEntries("events") << " events found" << endl;
  for (size_t i = 0; i < reader.getEntries("events"); i++) {// begin event loop
    const auto event = reader.readNextFrame("events");
    if (i%100==0) cout << i << " events processed" << endl;

    // Retrieve Inclusive Kinematics Collections
    auto& kin_truth = event.get<edm4eic::InclusiveKinematicsCollection>("InclusiveKinematicsTruth");
    auto& kin_electron = event.get<edm4eic::InclusiveKinematicsCollection>("InclusiveKinematicsElectron");
    auto& kin_jb = event.get<edm4eic::InclusiveKinematicsCollection>("InclusiveKinematicsJB");
    auto& kin_da = event.get<edm4eic::InclusiveKinematicsCollection>("InclusiveKinematicsDA");
    auto& kin_sigma = event.get<edm4eic::InclusiveKinematicsCollection>("InclusiveKinematicsSigma");
    auto& kin_esigma = event.get<edm4eic::InclusiveKinematicsCollection>("InclusiveKinematicsESigma");

    if (kin_truth.empty() || kin_electron.empty() || kin_jb.empty()) continue;

    x_truth = kin_truth.x()[0];
    x_electron = kin_electron.x()[0];
    x_jb = kin_jb.x()[0];
    x_da = kin_da.x()[0];
    x_sigma = kin_sigma.x()[0];
    x_esigma = kin_esigma.x()[0];

    y_truth = kin_truth.y()[0];
    y_electron = kin_electron.y()[0];
    y_jb = kin_jb.y()[0];
    y_da = kin_da.y()[0];
    y_sigma = kin_sigma.y()[0];
    y_esigma = kin_esigma.y()[0];

    Q2_truth = kin_truth.Q2()[0];
    Q2_electron = kin_electron.Q2()[0];
    Q2_jb = kin_jb.Q2()[0];
    Q2_da = kin_da.Q2()[0];
    Q2_sigma = kin_sigma.Q2()[0];
    Q2_esigma = kin_esigma.Q2()[0];

    // Some example cuts
    bool cuts = true;
    cuts = cuts && (y_truth < 0.95);
    cuts = cuts && (y_truth > 0.01);
    cuts = cuts && (Q2_truth > 1);

    if (!cuts) continue;

    hResoX_electron->Fill((x_electron-x_truth)/x_truth);
    hResoX_jb->Fill((x_jb-x_truth)/x_truth);
    hResoX_da->Fill((x_da-x_truth)/x_truth);
    hResoX_sigma->Fill((x_sigma-x_truth)/x_truth);
    hResoX_esigma->Fill((x_esigma-x_truth)/x_truth);

    hResoY_electron->Fill((y_electron-y_truth)/y_truth);
    hResoY_jb->Fill((y_jb-y_truth)/y_truth);
    hResoY_da->Fill((y_da-y_truth)/y_truth);
    hResoY_sigma->Fill((y_sigma-y_truth)/y_truth);
    hResoY_esigma->Fill((y_esigma-y_truth)/y_truth);

    hResoQ2_electron->Fill((Q2_electron-Q2_truth)/Q2_truth);
    hResoQ2_jb->Fill((Q2_jb-Q2_truth)/Q2_truth);
    hResoQ2_da->Fill((Q2_da-Q2_truth)/Q2_truth);
    hResoQ2_sigma->Fill((Q2_sigma-Q2_truth)/Q2_truth);
    hResoQ2_esigma->Fill((Q2_esigma-Q2_truth)/Q2_truth);

    hResoX_2D_electron->Fill(y_truth, (x_electron-x_truth)/x_truth);
    hResoX_2D_jb->Fill(y_truth, (x_jb-x_truth)/x_truth);
    hResoX_2D_da->Fill(y_truth, (x_da-x_truth)/x_truth);
    hResoX_2D_sigma->Fill(y_truth, (x_sigma-x_truth)/x_truth);
    hResoX_2D_esigma->Fill(y_truth, (x_esigma-x_truth)/x_truth);

    hResoY_2D_electron->Fill(y_truth, (y_electron-y_truth)/y_truth);
    hResoY_2D_jb->Fill(y_truth, (y_jb-y_truth)/y_truth);
    hResoY_2D_da->Fill(y_truth, (y_da-y_truth)/y_truth);
    hResoY_2D_sigma->Fill(y_truth, (y_sigma-y_truth)/y_truth);
    hResoY_2D_esigma->Fill(y_truth, (y_esigma-y_truth)/y_truth);

    hResoQ2_2D_electron->Fill(y_truth, (Q2_electron-Q2_truth)/Q2_truth);
    hResoQ2_2D_jb->Fill(y_truth, (Q2_jb-Q2_truth)/Q2_truth);
    hResoQ2_2D_da->Fill(y_truth, (Q2_da-Q2_truth)/Q2_truth);
    hResoQ2_2D_sigma->Fill(y_truth, (Q2_sigma-Q2_truth)/Q2_truth);
    hResoQ2_2D_esigma->Fill(y_truth, (Q2_esigma-Q2_truth)/Q2_truth);
  }// end event loop

  // Drawing the histograms
  auto canvas_x_1D = new TCanvas();
  canvas_x_1D->Divide(3,2);
  canvas_x_1D->cd(1);hResoX_electron->Draw("hist");
  canvas_x_1D->cd(2);hResoX_jb->Draw("hist");
  canvas_x_1D->cd(3);hResoX_da->Draw("hist");
  canvas_x_1D->cd(4);hResoX_sigma->Draw("hist");
  canvas_x_1D->cd(5);hResoX_esigma->Draw("hist");
  auto canvas_y_1D = new TCanvas();
  canvas_y_1D->Divide(3,2);
  canvas_y_1D->cd(1);hResoY_electron->Draw("hist");
  canvas_y_1D->cd(2);hResoY_jb->Draw("hist");
  canvas_y_1D->cd(3);hResoY_da->Draw("hist");
  canvas_y_1D->cd(4);hResoY_sigma->Draw("hist");
  canvas_y_1D->cd(5);hResoY_esigma->Draw("hist");
  auto canvas_Q2_1D = new TCanvas();
  canvas_Q2_1D->Divide(3,2);
  canvas_Q2_1D->cd(1);hResoQ2_electron->Draw("hist");
  canvas_Q2_1D->cd(2);hResoQ2_jb->Draw("hist");
  canvas_Q2_1D->cd(3);hResoQ2_da->Draw("hist");
  canvas_Q2_1D->cd(4);hResoQ2_sigma->Draw("hist");
  canvas_Q2_1D->cd(5);hResoQ2_esigma->Draw("hist");
  
  auto canvas_x_2D = new TCanvas();
  canvas_x_2D->Divide(3,2);
  canvas_x_2D->cd(1);if(bin_log) gPad->SetLogx();hResoX_2D_electron->Draw("colz");
  canvas_x_2D->cd(2);if(bin_log) gPad->SetLogx();hResoX_2D_jb->Draw("colz");
  canvas_x_2D->cd(3);if(bin_log) gPad->SetLogx();hResoX_2D_da->Draw("colz");
  canvas_x_2D->cd(4);if(bin_log) gPad->SetLogx();hResoX_2D_sigma->Draw("colz");
  canvas_x_2D->cd(5);if(bin_log) gPad->SetLogx();hResoX_2D_esigma->Draw("colz");
  auto canvas_y_2D = new TCanvas();
  canvas_y_2D->Divide(3,2);
  canvas_y_2D->cd(1);if(bin_log) gPad->SetLogx();hResoY_2D_electron->Draw("colz");
  canvas_y_2D->cd(2);if(bin_log) gPad->SetLogx();hResoY_2D_jb->Draw("colz");
  canvas_y_2D->cd(3);if(bin_log) gPad->SetLogx();hResoY_2D_da->Draw("colz");
  canvas_y_2D->cd(4);if(bin_log) gPad->SetLogx();hResoY_2D_sigma->Draw("colz");
  canvas_y_2D->cd(5);if(bin_log) gPad->SetLogx();hResoY_2D_esigma->Draw("colz");
  auto canvas_Q2_2D = new TCanvas();
  canvas_Q2_2D->Divide(3,2);
  canvas_Q2_2D->cd(1);if(bin_log) gPad->SetLogx();hResoQ2_2D_electron->Draw("colz");
  canvas_Q2_2D->cd(2);if(bin_log) gPad->SetLogx();hResoQ2_2D_jb->Draw("colz");
  canvas_Q2_2D->cd(3);if(bin_log) gPad->SetLogx();hResoQ2_2D_da->Draw("colz");
  canvas_Q2_2D->cd(4);if(bin_log) gPad->SetLogx();hResoQ2_2D_sigma->Draw("colz");
  canvas_Q2_2D->cd(5);if(bin_log) gPad->SetLogx();hResoQ2_2D_esigma->Draw("colz");

  cout << "Done!" << endl;
}

This script sets up the benchmark histograms, fills them in the event loop, and then draws them. Here, the resolutions on the reconstructed kinematic variables are chosen as the benchmarks, both the 1-dimensional (reco-true)/true distribution, and also 2-dimensional plots vs inelasticity, y. For a good reconstruction method, the (reco-true)/true distribution is centred on zero, with small fluctuations.

Run this script as

BASH

root -l BenchmarkReconstruction.C\(\"your_file.root\"\)

or as

BASH

root -l BenchmarkReconstruction.C\(\"your_file.root\",true\)

to bin logarithmically in inelasticity.

You may wish to investigate how the resolutions change in a scenario more relevant to your analysis. A set of example cuts are provided in the script

CPP

// Some example cuts
bool cuts = true;
cuts = cuts && (y_truth < 0.95);
cuts = cuts && (y_truth > 0.01);
cuts = cuts && (Q2_truth > 1);

These can be replaced with whatever cuts are used in your analysis, or you could use them to select areas of the phase space that you wish to investigate.

Challenge

Exercise

Modify the example cuts in BenchmarkReconstruction.C to isolate a region of phase space relevant to your analysis (for example a high-inelasticity or high-Q2 selection) and re-run the benchmark. Which reconstruction method gives the narrowest (reco-true)/true resolution in that region?

Show me the solution

There is no single correct answer - it depends on the region you select. In general the electron method performs best at high inelasticity y, while the JB and Double Angle methods do better at low y. The point of the exercise is to see that the “best” method is region-dependent, so you should benchmark the methods in the region relevant to your own analysis.

Key Points

• Use the podio Reader (podio::makeReader) to process simulation files using the data types implemented in edm4hep/edm4eic.

Content from Manual Reconstruction

---

Last updated on 2026-07-10 | Edit this page

Overview

Questions

• How do I reconstruct the inclusive kinematics myself?

Objectives

• Understand how to calculate the inclusive kinematics manually.
• Implement the various reconstruction methods in code.
• Verify your manual calculations through comparison to the InclusiveKinematicsXX values.

Doing the reconstruction yourself

---

It may be that you don’t want to use the default reconstruction provided in the InclusiveKinematicsXX branches - maybe you want to test a new electron finding or particle flow algorithm? In such cases you will need to calculate the kinematics yourself, as the InclusiveKinematicsXX branches only perform the reconstruction for a single scenario e.g. perfect electron ID, electron energy from tracking etc. If you’re having to write the reconstruction methods in your own code, it’s good to verify that your implementation of the methods is correct.

We can do this by comparing our manual calculations to the results stored in the InclusiveKinematicsXX branches. Copy the script below into a file called ManualReconstruction.C

CPP

// PODIO
#include "podio/Frame.h"
#include "podio/Reader.h"

// DATA MODEL
#include "edm4eic/InclusiveKinematicsCollection.h"
#include "edm4eic/ReconstructedParticleCollection.h"
#include "edm4eic/HadronicFinalStateCollection.h"
#include "edm4eic/ClusterCollection.h"
#include "edm4hep/Vector3f.h"

std::vector<float> calc_elec_method(float E, float theta, float pt_had, float sigma_h, float E_ebeam, float E_pbeam);
std::vector<float> calc_jb_method(float E, float theta, float pt_had, float sigma_h, float E_ebeam, float E_pbeam);
std::vector<float> calc_da_method(float E, float theta, float pt_had, float sigma_h, float E_ebeam, float E_pbeam);
std::vector<float> calc_sig_method(float E, float theta, float pt_had, float sigma_h, float E_ebeam, float E_pbeam);
std::vector<float> calc_esig_method(float E, float theta, float pt_had, float sigma_h, float E_ebeam, float E_pbeam);

void ManualReconstruction(std::string filename) {

  // Settings
  Float_t E_ebeam = 18;
  Float_t E_pbeam = 275;
  Float_t m_e = 0.000511;
  double xAngle = 25e-3;
  TLorentzVector pni, ei;
  ei.SetPxPyPzE(0, 0, -E_ebeam, E_ebeam);
  pni.SetPxPyPzE(-1*E_pbeam*TMath::Sin(xAngle), 0, E_pbeam*TMath::Cos(xAngle), E_pbeam);

  std::vector<std::string> inFiles = {filename};

  auto reader = podio::makeReader(inFiles);

  // Declare benchmark histograms
  TH1F *hResoX_electron = new TH1F("hResoX_electron","Electron method;#Deltax/x;Counts",500,-1,1);
  TH1F *hResoX_jb = new TH1F("hResoX_jb","JB method;#Deltax/x;Counts",500,-1,1);
  TH1F *hResoX_da = new TH1F("hResoX_da","Double Angle method;#Deltax/x;Counts",500,-1,1);
  TH1F *hResoX_sigma = new TH1F("hResoX_sigma","#Sigma method;#Deltax/x;Counts",500,-1,1);
  TH1F *hResoX_esigma = new TH1F("hResoX_esigma","e-#Sigma method;#Deltax/x;Counts",500,-1,1);

  TH1F *hResoY_electron = new TH1F("hResoY_electron","Electron method;#Deltay/y;Counts",500,-1,1);
  TH1F *hResoY_jb = new TH1F("hResoY_jb","JB method;#Deltay/y;Counts",500,-1,1);
  TH1F *hResoY_da = new TH1F("hResoY_da","Double Angle method;#Deltay/y;Counts",500,-1,1);
  TH1F *hResoY_sigma = new TH1F("hResoY_sigma","#Sigma method;#Deltay/y;Counts",500,-1,1);
  TH1F *hResoY_esigma = new TH1F("hResoY_esigma","e-#Sigma method;#Deltay/y;Counts",500,-1,1);

  TH1F *hResoQ2_electron = new TH1F("hResoQ2_electron","Electron method;#DeltaQ2/Q2;Counts",500,-1,1);
  TH1F *hResoQ2_jb = new TH1F("hResoQ2_jb","JB method;#DeltaQ2/Q2;Counts",500,-1,1);
  TH1F *hResoQ2_da = new TH1F("hResoQ2_da","Double Angle method;#DeltaQ2/Q2;Counts",500,-1,1);
  TH1F *hResoQ2_sigma = new TH1F("hResoQ2_sigma","#Sigma method;#DeltaQ2/Q2;Counts",500,-1,1);
  TH1F *hResoQ2_esigma = new TH1F("hResoQ2_esigma","e-#Sigma method;#DeltaQ2/Q2;Counts",500,-1,1);

  Float_t x_truth, x_electron, x_jb, x_da, x_sigma, x_esigma;
  Float_t y_truth, y_electron, y_jb, y_da, y_sigma, y_esigma;
  Float_t Q2_truth, Q2_electron, Q2_jb, Q2_da, Q2_sigma, Q2_esigma;
  Float_t E, theta, sigma_h, pt_had;

  cout << reader.getEntries("events") << " events found" << endl;
  for (size_t i = 0; i < reader.getEntries("events"); i++) {// begin event loop
    const auto event = reader.readNextFrame("events");
    if (i%100==0) cout << i << " events processed" << endl;

    // Retrieve Inclusive Kinematics Collections
    auto& kin_truth = event.get<edm4eic::InclusiveKinematicsCollection>("InclusiveKinematicsTruth");
    auto& kin_electron = event.get<edm4eic::InclusiveKinematicsCollection>("InclusiveKinematicsElectron");
    auto& kin_jb = event.get<edm4eic::InclusiveKinematicsCollection>("InclusiveKinematicsJB");
    auto& kin_da = event.get<edm4eic::InclusiveKinematicsCollection>("InclusiveKinematicsDA");
    auto& kin_sigma = event.get<edm4eic::InclusiveKinematicsCollection>("InclusiveKinematicsSigma");
    auto& kin_esigma = event.get<edm4eic::InclusiveKinematicsCollection>("InclusiveKinematicsESigma");

    // Retrieve Scattered electron and HFS
    auto& eleCollection = event.get<edm4eic::ReconstructedParticleCollection>("ScatteredElectronsTruth");
    auto& hfsCollection = event.get<edm4eic::HadronicFinalStateCollection>("HadronicFinalState");

    // Store kinematics from InclusiveKinematics branches
    if (kin_truth.empty() || kin_electron.empty() || kin_jb.empty()) continue;
    x_truth = kin_truth.x()[0];
    x_electron = kin_electron.x()[0];
    x_jb = kin_jb.x()[0];
    x_da = kin_da.x()[0];
    x_sigma = kin_sigma.x()[0];
    x_esigma = kin_esigma.x()[0];

    y_truth = kin_truth.y()[0];
    y_electron = kin_electron.y()[0];
    y_jb = kin_jb.y()[0];
    y_da = kin_da.y()[0];
    y_sigma = kin_sigma.y()[0];
    y_esigma = kin_esigma.y()[0];

    Q2_truth = kin_truth.Q2()[0];
    Q2_electron = kin_electron.Q2()[0];
    Q2_jb = kin_jb.Q2()[0];
    Q2_da = kin_da.Q2()[0];
    Q2_sigma = kin_sigma.Q2()[0];
    Q2_esigma = kin_esigma.Q2()[0];

    TLorentzVector scat_ele;
    E = eleCollection[0].getEnergy();
    auto& ele_momentum = eleCollection[0].getMomentum();
    scat_ele.SetPxPyPzE(ele_momentum.x, ele_momentum.y, ele_momentum.z, E);
    theta = scat_ele.Theta();
    sigma_h = hfsCollection[0].getSigma();
    pt_had = hfsCollection[0].getPT();

    // Calculate kinematics manually
    std::vector<float> elec_reco = calc_elec_method(E, theta, pt_had, sigma_h, E_ebeam, E_pbeam);
    std::vector<float> jb_reco = calc_jb_method(E, theta, pt_had, sigma_h, E_ebeam, E_pbeam);
    std::vector<float> da_reco = calc_da_method(E, theta, pt_had, sigma_h, E_ebeam, E_pbeam);
    std::vector<float> sigma_reco = calc_sig_method(E, theta, pt_had, sigma_h, E_ebeam, E_pbeam);
    std::vector<float> esigma_reco = calc_esig_method(E, theta, pt_had, sigma_h, E_ebeam, E_pbeam);

    // Some example cuts
    bool cuts = true;
    cuts = cuts && (y_truth < 0.95);
    cuts = cuts && (y_truth > 0.01);
    cuts = cuts && (Q2_truth > 1);

    if (!cuts) continue;

    // Fill histograms with difference of calculated kinematics
    // and those retrieved from InclusiveKinematics branches
    hResoX_electron->Fill((x_electron-elec_reco[0])/elec_reco[0]);
    hResoX_jb->Fill((x_jb-jb_reco[0])/jb_reco[0]);
    hResoX_da->Fill((x_da-da_reco[0])/da_reco[0]);
    hResoX_sigma->Fill((x_sigma-sigma_reco[0])/sigma_reco[0]);
    hResoX_esigma->Fill((x_esigma-esigma_reco[0])/esigma_reco[0]);

    hResoY_electron->Fill((y_electron-elec_reco[1])/elec_reco[1]);
    hResoY_jb->Fill((y_jb-jb_reco[1])/jb_reco[1]);
    hResoY_da->Fill((y_da-da_reco[1])/da_reco[1]);
    hResoY_sigma->Fill((y_sigma-sigma_reco[1])/sigma_reco[1]);
    hResoY_esigma->Fill((y_esigma-esigma_reco[1])/esigma_reco[1]);

    hResoQ2_electron->Fill((Q2_electron-elec_reco[2])/elec_reco[2]);
    hResoQ2_jb->Fill((Q2_jb-jb_reco[2])/jb_reco[2]);
    hResoQ2_da->Fill((Q2_da-da_reco[2])/da_reco[2]);
    hResoQ2_sigma->Fill((Q2_sigma-sigma_reco[2])/sigma_reco[2]);
    hResoQ2_esigma->Fill((Q2_esigma-esigma_reco[2])/esigma_reco[2]);
  }
  // Drawing the histograms
  auto canvas_x_1D = new TCanvas();
  canvas_x_1D->Divide(3,2);
  canvas_x_1D->cd(1);hResoX_electron->Draw("hist");
  canvas_x_1D->cd(2);hResoX_jb->Draw("hist");
  canvas_x_1D->cd(3);hResoX_da->Draw("hist");
  canvas_x_1D->cd(4);hResoX_sigma->Draw("hist");
  canvas_x_1D->cd(5);hResoX_esigma->Draw("hist");
  auto canvas_y_1D = new TCanvas();
  canvas_y_1D->Divide(3,2);
  canvas_y_1D->cd(1);hResoY_electron->Draw("hist");
  canvas_y_1D->cd(2);hResoY_jb->Draw("hist");
  canvas_y_1D->cd(3);hResoY_da->Draw("hist");
  canvas_y_1D->cd(4);hResoY_sigma->Draw("hist");
  canvas_y_1D->cd(5);hResoY_esigma->Draw("hist");
  auto canvas_Q2_1D = new TCanvas();
  canvas_Q2_1D->Divide(3,2);
  canvas_Q2_1D->cd(1);hResoQ2_electron->Draw("hist");
  canvas_Q2_1D->cd(2);hResoQ2_jb->Draw("hist");
  canvas_Q2_1D->cd(3);hResoQ2_da->Draw("hist");
  canvas_Q2_1D->cd(4);hResoQ2_sigma->Draw("hist");
  canvas_Q2_1D->cd(5);hResoQ2_esigma->Draw("hist");

  cout << "Done!" << endl;
}

// electron method
std::vector<float> calc_elec_method(float E, float theta, float pt_had, float sigma_h, float E_ebeam, float E_pbeam) {
  float Q2  = 2.*E_ebeam*E*(1+TMath::Cos(theta));
  float y = 1. - (E/E_ebeam)*TMath::Sin(theta/2)*TMath::Sin(theta/2);
  float x = Q2/(4*E_ebeam*E_pbeam*y);
  return {x, y, Q2};
}

// jb method
std::vector<float> calc_jb_method(float E, float theta, float pt_had, float sigma_h, float E_ebeam, float E_pbeam) {
  float y = sigma_h/(2*E_ebeam);
  float Q2 = pt_had*pt_had / (1-y);
  float x = Q2/(4*E_ebeam*E_pbeam*y);
  return {x, y, Q2};
}

// float angle method
std::vector<float> calc_da_method(float E, float theta, float pt_had, float sigma_h, float E_ebeam, float E_pbeam) {
  float alpha_h = sigma_h/pt_had;
  float alpha_e = TMath::Tan(theta/2);
  float y = alpha_h / (alpha_e + alpha_h);
  float Q2 = 4*E_ebeam*E_ebeam / (alpha_e * (alpha_h + alpha_e));
  float x = Q2/(4*E_ebeam*E_pbeam*y);
  return {x, y, Q2};
}

// sigma method
std::vector<float> calc_sig_method(float E, float theta, float pt_had, float sigma_h, float E_ebeam, float E_pbeam) {
  float y = sigma_h/(sigma_h + E*(1 - TMath::Cos(theta))); 
  float Q2 = E*E*TMath::Sin(theta)*TMath::Sin(theta) / (1-y);
  float x = Q2/(4*E_ebeam*E_pbeam*y);
  return {x, y, Q2};
}

// e-sigma method
std::vector<float> calc_esig_method(float E, float theta, float pt_had, float sigma_h, float E_ebeam, float E_pbeam) {
  float Q2  = 2.*E_ebeam*E*(1+TMath::Cos(theta));
  float x = calc_sig_method(E,theta,pt_had,sigma_h,E_ebeam,E_pbeam)[0];
  float y = Q2/(4*E_ebeam*E_pbeam*x);
  return {x, y, Q2};
}

As previously, you can run this script as

BASH

root -l ManualReconstruction.C\(\"your_file.root\"\)

This produces plots comparing the manual calculations to the values in the branches as (branch_calc-manual_calc)/manual_calc. The manual calculations were coded as

CPP

// electron method
std::vector<float> calc_elec_method(float E, float theta, float pt_had, float sigma_h, float E_ebeam, float E_pbeam) {
  float Q2  = 2.*E_ebeam*E*(1+TMath::Cos(theta));
  float y = 1. - (E/E_ebeam)*TMath::Sin(theta/2)*TMath::Sin(theta/2);
  float x = Q2/(4*E_ebeam*E_pbeam*y);
  return {x, y, Q2};
}

// jb method
std::vector<float> calc_jb_method(float E, float theta, float pt_had, float sigma_h, float E_ebeam, float E_pbeam) {
  float y = sigma_h/(2*E_ebeam);
  float Q2 = pt_had*pt_had / (1-y);
  float x = Q2/(4*E_ebeam*E_pbeam*y);
  return {x, y, Q2};
}

// float angle method
std::vector<float> calc_da_method(float E, float theta, float pt_had, float sigma_h, float E_ebeam, float E_pbeam) {
  float alpha_h = sigma_h/pt_had;
  float alpha_e = TMath::Tan(theta/2);
  float y = alpha_h / (alpha_e + alpha_h);
  float Q2 = 4*E_ebeam*E_ebeam / (alpha_e * (alpha_h + alpha_e));
  float x = Q2/(4*E_ebeam*E_pbeam*y);
  return {x, y, Q2};
}

// sigma method
std::vector<float> calc_sig_method(float E, float theta, float pt_had, float sigma_h, float E_ebeam, float E_pbeam) {
  float y = sigma_h/(sigma_h + E*(1 - TMath::Cos(theta))); 
  float Q2 = E*E*TMath::Sin(theta)*TMath::Sin(theta) / (1-y);
  float x = Q2/(4*E_ebeam*E_pbeam*y);
  return {x, y, Q2};
}

// e-sigma method
std::vector<float> calc_esig_method(float E, float theta, float pt_had, float sigma_h, float E_ebeam, float E_pbeam) {
  float Q2  = 2.*E_ebeam*E*(1+TMath::Cos(theta));
  float x = calc_sig_method(E,theta,pt_had,sigma_h,E_ebeam,E_pbeam)[0];
  float y = Q2/(4*E_ebeam*E_pbeam*x);
  return {x, y, Q2};
}

such that they take the basic quantities: the scattered electron energy and angle, the Hadronic Final State transverse momentum and E-pz sum, and the energies of the beam electrons/protons as inputs. We would expect these to give the exact same result as those produced by the InclusiveKinematicsXX branches, assuming that they are implemented the same way. For files in the March 2025 campaign, we see this to be case for all methods - except for the x and y calculations in the electron method. To understand this better, we can compare what is implemented in our version of the method to the calculations in the InclusiveKinematicsElectron algorithm.

Looking at this, we see that the InclusiveKinematicsElectron branch uses four vectors in its calculations, while the manual calculation here uses only the electron energy and angle, and the beam energies. Ordinarily, we would expect these to give equivalent results, if not for the fact that at ePIC there is a 25mrad crossing angle: we may therefore see different results for the Lorentz Invariant four vector based approach compared to the manual calculation, the equations for which are derived assuming head-on collisions.

Challenge

Exercise

Try implementing your own version of the electron method using four vectors, accounting for the 25mrad crossing angle. Verify that it matches the output of the InclusiveKinematicsElectron branch.

Show me the solution

Build the scattered-electron and beam four vectors (as pni and ei are set up at the top of ManualReconstruction.C, using the crossing angle) and compute Q2, y, and x from the Lorentz invariants (Q2 = -(ei - scat_ele)^2, y = (pni . (ei - scat_ele)) / (pni . ei), x = Q2 / (2 pni . (ei - scat_ele))). Because this treatment includes the crossing angle, the (branch - manual)/manual distributions for x and y should now be centred on zero, matching the InclusiveKinematicsElectron branch, unlike the head-on formula used above.

Key Points

• The reconstruction methods all use some combination of the same basic information: Ee, theta_e, sigma_h, pt_h.

Content from Optimise Reconstruction

---

Last updated on 2026-07-10 | Edit this page

Overview

Questions

• How can I improve the reconstruction when using standard methods?

Objectives

• Compare kinematic resolutions for electron reconstruction based on tracks alone and for tracks+calorimetry.
• Use realistic scattered electron ID in the reconstruction.

Optimising the reconstruction

---

There are four quantities that are used to reconstruct the inclusive kinematics: the scattered electron energy and polar angle, the hadronic final state (HFS) transverse momentum, and the E-pz sum of all particles in the HFS. The values of these quantities depends on the information used in the reconstruction of the scattered electron and the HFS. In some regions of the phase space, the calorimeters may do a much better job of reconstructing the scattered electron compared to the tracking system.

We can investigate this using the script below, which should be copied into a file called OptimiseReconstruction.C

CPP

// PODIO
#include "podio/Frame.h"
#include "podio/Reader.h"

// DATA MODEL
#include "edm4eic/InclusiveKinematicsCollection.h"
#include "edm4hep/MCParticleCollection.h"
#include "edm4eic/ReconstructedParticleCollection.h"
#include "edm4eic/MCRecoClusterParticleAssociationCollection.h"
#include "edm4eic/MCRecoParticleAssociationCollection.h"
#include "edm4eic/HadronicFinalStateCollection.h"
#include "edm4eic/ClusterCollection.h"
#include "edm4hep/Vector3f.h"

std::vector<float> calc_elec_method(float E, float theta, float pt_had, float sigma_h, float E_ebeam, float E_pbeam);
std::vector<float> calc_jb_method(float E, float theta, float pt_had, float sigma_h, float E_ebeam, float E_pbeam);
std::vector<float> calc_da_method(float E, float theta, float pt_had, float sigma_h, float E_ebeam, float E_pbeam);
std::vector<float> calc_sig_method(float E, float theta, float pt_had, float sigma_h, float E_ebeam, float E_pbeam);
std::vector<float> calc_esig_method(float E, float theta, float pt_had, float sigma_h, float E_ebeam, float E_pbeam);

void OptimiseReconstruction(std::string filename) {

  // Settings
  Float_t E_ebeam = 18;
  Float_t E_pbeam = 275;
  Float_t m_e = 0.000511;

  std::vector<std::string> inFiles = {filename};

  auto reader = podio::makeReader(inFiles);

  // Declare benchmark histograms
  TH1F *hResoX_electron = new TH1F("hResoX_electron","Electron method;#Deltax/x;Counts",100,-1,1);
  TH1F *hResoX_jb = new TH1F("hResoX_jb","JB method;#Deltax/x;Counts",100,-1,1);
  TH1F *hResoX_da = new TH1F("hResoX_da","Double Angle method;#Deltax/x;Counts",100,-1,1);
  TH1F *hResoX_sigma = new TH1F("hResoX_sigma","#Sigma method;#Deltax/x;Counts",100,-1,1);
  TH1F *hResoX_esigma = new TH1F("hResoX_esigma","e-#Sigma method;#Deltax/x;Counts",100,-1,1);

  TH1F *hResoY_electron = new TH1F("hResoY_electron","Electron method;#Deltay/y;Counts",100,-1,1);
  TH1F *hResoY_jb = new TH1F("hResoY_jb","JB method;#Deltay/y;Counts",100,-1,1);
  TH1F *hResoY_da = new TH1F("hResoY_da","Double Angle method;#Deltay/y;Counts",100,-1,1);
  TH1F *hResoY_sigma = new TH1F("hResoY_sigma","#Sigma method;#Deltay/y;Counts",100,-1,1);
  TH1F *hResoY_esigma = new TH1F("hResoY_esigma","e-#Sigma method;#Deltay/y;Counts",100,-1,1);

  TH1F *hResoQ2_electron = new TH1F("hResoQ2_electron","Electron method;#DeltaQ2/Q2;Counts",100,-1,1);
  TH1F *hResoQ2_jb = new TH1F("hResoQ2_jb","JB method;#DeltaQ2/Q2;Counts",100,-1,1);
  TH1F *hResoQ2_da = new TH1F("hResoQ2_da","Double Angle method;#DeltaQ2/Q2;Counts",100,-1,1);
  TH1F *hResoQ2_sigma = new TH1F("hResoQ2_sigma","#Sigma method;#DeltaQ2/Q2;Counts",100,-1,1);
  TH1F *hResoQ2_esigma = new TH1F("hResoQ2_esigma","e-#Sigma method;#DeltaQ2/Q2;Counts",100,-1,1);

  TH1F *hResoX_calo_electron = new TH1F("hResoX_calo_electron","Electron method;#Deltax/x;Counts",100,-1,1);
  TH1F *hResoX_calo_jb = new TH1F("hResoX_calo_jb","JB method;#Deltax/x;Counts",100,-1,1);
  TH1F *hResoX_calo_da = new TH1F("hResoX_calo_da","Double Angle method;#Deltax/x;Counts",100,-1,1);
  TH1F *hResoX_calo_sigma = new TH1F("hResoX_calo_sigma","#Sigma method;#Deltax/x;Counts",100,-1,1);
  TH1F *hResoX_calo_esigma = new TH1F("hResoX_calo_esigma","e-#Sigma method;#Deltax/x;Counts",100,-1,1);

  TH1F *hResoY_calo_electron = new TH1F("hResoY_calo_electron","Electron method;#Deltay/y;Counts",100,-1,1);
  TH1F *hResoY_calo_jb = new TH1F("hResoY_calo_jb","JB method;#Deltay/y;Counts",100,-1,1);
  TH1F *hResoY_calo_da = new TH1F("hResoY_calo_da","Double Angle method;#Deltay/y;Counts",100,-1,1);
  TH1F *hResoY_calo_sigma = new TH1F("hResoY_calo_sigma","#Sigma method;#Deltay/y;Counts",100,-1,1);
  TH1F *hResoY_calo_esigma = new TH1F("hResoY_calo_esigma","e-#Sigma method;#Deltay/y;Counts",100,-1,1);

  TH1F *hResoQ2_calo_electron = new TH1F("hResoQ2_calo_electron","Electron method;#DeltaQ2/Q2;Counts",100,-1,1);
  TH1F *hResoQ2_calo_jb = new TH1F("hResoQ2_calo_jb","JB method;#DeltaQ2/Q2;Counts",100,-1,1);
  TH1F *hResoQ2_calo_da = new TH1F("hResoQ2_calo_da","Double Angle method;#DeltaQ2/Q2;Counts",100,-1,1);
  TH1F *hResoQ2_calo_sigma = new TH1F("hResoQ2_calo_sigma","#Sigma method;#DeltaQ2/Q2;Counts",100,-1,1);
  TH1F *hResoQ2_calo_esigma = new TH1F("hResoQ2_calo_esigma","e-#Sigma method;#DeltaQ2/Q2;Counts",100,-1,1);

  Float_t x_truth, x_electron, x_jb, x_da, x_sigma, x_esigma;
  Float_t y_truth, y_electron, y_jb, y_da, y_sigma, y_esigma;
  Float_t Q2_truth, Q2_electron, Q2_jb, Q2_da, Q2_sigma, Q2_esigma;
  Float_t E, theta, sigma_h, pt_had;

  cout << reader.getEntries("events") << " events found" << endl;
  for (size_t i = 0; i < reader.getEntries("events"); i++) {// begin event loop
    const auto event = reader.readNextFrame("events");
    if (i%100==0) cout << i << " events processed" << endl;

    // Retrieve Inclusive Kinematics Collections
    auto& kin_truth = event.get<edm4eic::InclusiveKinematicsCollection>("InclusiveKinematicsTruth");
    auto& kin_electron = event.get<edm4eic::InclusiveKinematicsCollection>("InclusiveKinematicsElectron");
    auto& kin_jb = event.get<edm4eic::InclusiveKinematicsCollection>("InclusiveKinematicsJB");
    auto& kin_da = event.get<edm4eic::InclusiveKinematicsCollection>("InclusiveKinematicsDA");
    auto& kin_sigma = event.get<edm4eic::InclusiveKinematicsCollection>("InclusiveKinematicsSigma");
    auto& kin_esigma = event.get<edm4eic::InclusiveKinematicsCollection>("InclusiveKinematicsESigma");

    // Retrieve Scattered electron and HFS
    auto& eleCollection = event.get<edm4eic::ReconstructedParticleCollection>("ScatteredElectronsEMinusPz");
    auto& hfsCollection = event.get<edm4eic::HadronicFinalStateCollection>("HadronicFinalState");

    // Store kinematics from InclusiveKinematics branches
    if (kin_truth.empty() || kin_electron.empty() || kin_jb.empty()) continue;
    if (eleCollection.empty()) continue;

    TLorentzVector scat_ele;
    E = eleCollection[0].getEnergy();
    auto& ele_momentum = eleCollection[0].getMomentum();
    scat_ele.SetPxPyPzE(ele_momentum.x, ele_momentum.y, ele_momentum.z, E);
    theta = scat_ele.Theta();
    sigma_h = hfsCollection[0].getSigma();
    pt_had = hfsCollection[0].getPT();

    // Calculate kinematics manually
    std::vector<float> elec_reco = calc_elec_method(E, theta, pt_had, sigma_h, E_ebeam, E_pbeam);
    std::vector<float> jb_reco = calc_jb_method(E, theta, pt_had, sigma_h, E_ebeam, E_pbeam);
    std::vector<float> da_reco = calc_da_method(E, theta, pt_had, sigma_h, E_ebeam, E_pbeam);
    std::vector<float> sigma_reco = calc_sig_method(E, theta, pt_had, sigma_h, E_ebeam, E_pbeam);
    std::vector<float> esigma_reco = calc_esig_method(E, theta, pt_had, sigma_h, E_ebeam, E_pbeam);

    x_truth = kin_truth.x()[0];
    y_truth = kin_truth.y()[0];
    Q2_truth = kin_truth.Q2()[0];

    x_electron = elec_reco[0];
    x_jb = jb_reco[0];
    x_da = da_reco[0];
    x_sigma = sigma_reco[0];
    x_esigma = esigma_reco[0];

    y_electron = elec_reco[1];
    y_jb = jb_reco[1];
    y_da = da_reco[1];
    y_sigma = sigma_reco[1];
    y_esigma = esigma_reco[1];

    Q2_electron = elec_reco[2];
    Q2_jb = jb_reco[2];
    Q2_da = da_reco[2];
    Q2_sigma = sigma_reco[2];
    Q2_esigma = esigma_reco[2];

    // Some example cuts
    bool cuts = true;
    cuts = cuts && (y_truth < 0.95);
    cuts = cuts && (y_truth > 0.01);
    cuts = cuts && (Q2_truth > 1);

    if (!cuts) continue;

    hResoX_electron->Fill((x_electron-x_truth)/x_truth);
    hResoX_jb->Fill((x_jb-x_truth)/x_truth);
    hResoX_da->Fill((x_da-x_truth)/x_truth);
    hResoX_sigma->Fill((x_sigma-x_truth)/x_truth);
    hResoX_esigma->Fill((x_esigma-x_truth)/x_truth);

    hResoY_electron->Fill((y_electron-y_truth)/y_truth);
    hResoY_jb->Fill((y_jb-y_truth)/y_truth);
    hResoY_da->Fill((y_da-y_truth)/y_truth);
    hResoY_sigma->Fill((y_sigma-y_truth)/y_truth);
    hResoY_esigma->Fill((y_esigma-y_truth)/y_truth);

    hResoQ2_electron->Fill((Q2_electron-Q2_truth)/Q2_truth);
    hResoQ2_jb->Fill((Q2_jb-Q2_truth)/Q2_truth);
    hResoQ2_da->Fill((Q2_da-Q2_truth)/Q2_truth);
    hResoQ2_sigma->Fill((Q2_sigma-Q2_truth)/Q2_truth);
    hResoQ2_esigma->Fill((Q2_esigma-Q2_truth)/Q2_truth);

    // Replace electron momentum with energy from calo cluster
    E = eleCollection[0].getClusters()[0].getEnergy();

    // Recalculate kinematics
    std::vector<float> calo_elec_reco = calc_elec_method(E, theta, pt_had, sigma_h, E_ebeam, E_pbeam);
    std::vector<float> calo_jb_reco = calc_jb_method(E, theta, pt_had, sigma_h, E_ebeam, E_pbeam);
    std::vector<float> calo_da_reco = calc_da_method(E, theta, pt_had, sigma_h, E_ebeam, E_pbeam);
    std::vector<float> calo_sigma_reco = calc_sig_method(E, theta, pt_had, sigma_h, E_ebeam, E_pbeam);
    std::vector<float> calo_esigma_reco = calc_esig_method(E, theta, pt_had, sigma_h, E_ebeam, E_pbeam);

    x_electron = calo_elec_reco[0];
    x_jb = calo_jb_reco[0];
    x_da = calo_da_reco[0];
    x_sigma = calo_sigma_reco[0];
    x_esigma = calo_esigma_reco[0];

    y_electron = calo_elec_reco[1];
    y_jb = calo_jb_reco[1];
    y_da = calo_da_reco[1];
    y_sigma = calo_sigma_reco[1];
    y_esigma = calo_esigma_reco[1];

    Q2_electron = calo_elec_reco[2];
    Q2_jb = calo_jb_reco[2];
    Q2_da = calo_da_reco[2];
    Q2_sigma = calo_sigma_reco[2];
    Q2_esigma = calo_esigma_reco[2];

    hResoX_calo_electron->Fill((x_electron-x_truth)/x_truth);
    hResoX_calo_jb->Fill((x_jb-x_truth)/x_truth);
    hResoX_calo_da->Fill((x_da-x_truth)/x_truth);
    hResoX_calo_sigma->Fill((x_sigma-x_truth)/x_truth);
    hResoX_calo_esigma->Fill((x_esigma-x_truth)/x_truth);

    hResoY_calo_electron->Fill((y_electron-y_truth)/y_truth);
    hResoY_calo_jb->Fill((y_jb-y_truth)/y_truth);
    hResoY_calo_da->Fill((y_da-y_truth)/y_truth);
    hResoY_calo_sigma->Fill((y_sigma-y_truth)/y_truth);
    hResoY_calo_esigma->Fill((y_esigma-y_truth)/y_truth);

    hResoQ2_calo_electron->Fill((Q2_electron-Q2_truth)/Q2_truth);
    hResoQ2_calo_jb->Fill((Q2_jb-Q2_truth)/Q2_truth);
    hResoQ2_calo_da->Fill((Q2_da-Q2_truth)/Q2_truth);
    hResoQ2_calo_sigma->Fill((Q2_sigma-Q2_truth)/Q2_truth);
    hResoQ2_calo_esigma->Fill((Q2_esigma-Q2_truth)/Q2_truth);
  }
  // Drawing the histograms
  auto canvas_x_1D = new TCanvas();
  canvas_x_1D->Divide(3,2);
  canvas_x_1D->cd(1); hResoX_electron->Draw("hist");
  hResoX_calo_electron->SetLineStyle(2); hResoX_calo_electron->SetLineColor(kRed); hResoX_calo_electron->Draw("hist same");
  canvas_x_1D->cd(2); hResoX_jb->Draw("hist");
  hResoX_calo_jb->SetLineStyle(2); hResoX_calo_jb->SetLineColor(kRed); hResoX_calo_jb->Draw("hist same");
  canvas_x_1D->cd(3); hResoX_da->Draw("hist");
  hResoX_calo_da->SetLineStyle(2); hResoX_calo_da->SetLineColor(kRed); hResoX_calo_da->Draw("hist same");
  canvas_x_1D->cd(4); hResoX_sigma->Draw("hist");
  hResoX_calo_sigma->SetLineStyle(2); hResoX_calo_sigma->SetLineColor(kRed); hResoX_calo_sigma->Draw("hist same");
  canvas_x_1D->cd(5); hResoX_esigma->Draw("hist");
  hResoX_calo_esigma->SetLineStyle(2); hResoX_calo_esigma->SetLineColor(kRed); hResoX_calo_esigma->Draw("hist same");
  canvas_x_1D->cd(6);
  auto legend_x = new TLegend(0.1,0.1,0.9,0.9);
  legend_x->AddEntry(hResoX_electron, "E_{e} from tracker", "l");
  legend_x->AddEntry(hResoX_calo_electron, "E_{e} from ECAL", "l");
  legend_x->Draw();
  
  auto canvas_y_1D = new TCanvas();
  canvas_y_1D->Divide(3,2);
  canvas_y_1D->cd(1); hResoY_electron->Draw("hist");
  hResoY_calo_electron->SetLineStyle(2); hResoY_calo_electron->SetLineColor(kRed); hResoY_calo_electron->Draw("hist same");
  canvas_y_1D->cd(2); hResoY_jb->Draw("hist");
  hResoY_calo_jb->SetLineStyle(2); hResoY_calo_jb->SetLineColor(kRed); hResoY_calo_jb->Draw("hist same");
  canvas_y_1D->cd(3); hResoY_da->Draw("hist");
  hResoY_calo_da->SetLineStyle(2); hResoY_calo_da->SetLineColor(kRed); hResoY_calo_da->Draw("hist same");
  canvas_y_1D->cd(4); hResoY_sigma->Draw("hist");
  hResoY_calo_sigma->SetLineStyle(2); hResoY_calo_sigma->SetLineColor(kRed); hResoY_calo_sigma->Draw("hist same");
  canvas_y_1D->cd(5); hResoY_esigma->Draw("hist");
  hResoY_calo_esigma->SetLineStyle(2); hResoY_calo_esigma->SetLineColor(kRed); hResoY_calo_esigma->Draw("hist same");
  canvas_y_1D->cd(6);
  auto legend_y = new TLegend(0.1,0.1,0.9,0.9);
  legend_y->AddEntry(hResoY_electron, "E_{e} from tracker", "l");
  legend_y->AddEntry(hResoY_calo_electron, "E_{e} from ECAL", "l");
  legend_y->Draw();
  
  auto canvas_Q2_1D = new TCanvas();
  canvas_Q2_1D->Divide(3,2);
  canvas_Q2_1D->cd(1); hResoQ2_electron->Draw("hist");
  hResoQ2_calo_electron->SetLineStyle(2); hResoQ2_calo_electron->SetLineColor(kRed); hResoQ2_calo_electron->Draw("hist same");
  canvas_Q2_1D->cd(2); hResoQ2_jb->Draw("hist");
  hResoQ2_calo_jb->SetLineStyle(2); hResoQ2_calo_jb->SetLineColor(kRed); hResoQ2_calo_jb->Draw("hist same");
  canvas_Q2_1D->cd(3); hResoQ2_da->Draw("hist");
  hResoQ2_calo_da->SetLineStyle(2); hResoQ2_calo_da->SetLineColor(kRed); hResoQ2_calo_da->Draw("hist same");
  canvas_Q2_1D->cd(4); hResoQ2_sigma->Draw("hist");
  hResoQ2_calo_sigma->SetLineStyle(2); hResoQ2_calo_sigma->SetLineColor(kRed); hResoQ2_calo_sigma->Draw("hist same");
  canvas_Q2_1D->cd(5); hResoQ2_esigma->Draw("hist");
  hResoQ2_calo_esigma->SetLineStyle(2); hResoQ2_calo_esigma->SetLineColor(kRed); hResoQ2_calo_esigma->Draw("hist same");
  canvas_Q2_1D->cd(6);
  auto legend_Q2 = new TLegend(0.1,0.1,0.9,0.9);
  legend_Q2->AddEntry(hResoQ2_electron, "E_{e} from tracker", "l");
  legend_Q2->AddEntry(hResoQ2_calo_electron, "E_{e} from ECAL", "l");
  legend_Q2->Draw();
  

  cout << "Done!" << endl;
}

// electron method
std::vector<float> calc_elec_method(float E, float theta, float pt_had, float sigma_h, float E_ebeam, float E_pbeam) {
  float Q2  = 2.*E_ebeam*E*(1+TMath::Cos(theta));
  float y = 1. - (E/E_ebeam)*TMath::Sin(theta/2)*TMath::Sin(theta/2);
  float x = Q2/(4*E_ebeam*E_pbeam*y);
  return {x, y, Q2};
}

// jb method
std::vector<float> calc_jb_method(float E, float theta, float pt_had, float sigma_h, float E_ebeam, float E_pbeam) {
  float y = sigma_h/(2*E_ebeam);
  float Q2 = pt_had*pt_had / (1-y);
  float x = Q2/(4*E_ebeam*E_pbeam*y);
  return {x, y, Q2};
}

// float angle method
std::vector<float> calc_da_method(float E, float theta, float pt_had, float sigma_h, float E_ebeam, float E_pbeam) {
  float alpha_h = sigma_h/pt_had;
  float alpha_e = TMath::Tan(theta/2);
  float y = alpha_h / (alpha_e + alpha_h);
  float Q2 = 4*E_ebeam*E_ebeam / (alpha_e * (alpha_h + alpha_e));
  float x = Q2/(4*E_ebeam*E_pbeam*y);
  return {x, y, Q2};
}

// sigma method
std::vector<float> calc_sig_method(float E, float theta, float pt_had, float sigma_h, float E_ebeam, float E_pbeam) {
  float y = sigma_h/(sigma_h + E*(1 - TMath::Cos(theta))); 
  float Q2 = E*E*TMath::Sin(theta)*TMath::Sin(theta) / (1-y);
  float x = Q2/(4*E_ebeam*E_pbeam*y);
  return {x, y, Q2};
}

// e-sigma method
std::vector<float> calc_esig_method(float E, float theta, float pt_had, float sigma_h, float E_ebeam, float E_pbeam) {
  float Q2  = 2.*E_ebeam*E*(1+TMath::Cos(theta));
  float x = calc_sig_method(E,theta,pt_had,sigma_h,E_ebeam,E_pbeam)[0];
  float y = Q2/(4*E_ebeam*E_pbeam*x);
  return {x, y, Q2};
}

This script combines features of the scripts shown in the previous two sections. The plots that are produced when you run

BASH

root -l OptimiseReconstruction.C\(\"your_file.root\"\)

show the (reco-true)/true distributions as before, with the reconstructed values coming from the manual calculations.

The output of the basic electron-finder implemented in EICrecon is found in the ScatteredElectronsEMinusPz branch, which is accessed as

CPP

auto& eleCollection = event.get<edm4eic::ReconstructedParticleCollection>("ScatteredElectronsEMinusPz");

This returns a list of all particles in ReconstructedParticles that have matched tracks and ECAL clusters that pass an E/p cut, ordered by momentum. In this code, we take the first element (largest momentum) of the list - be aware that this assumption causes a drop in efficiency larger values of inelasticity.

In this script we compare the quality of the reconstruction methods for two different scenarios: the first where the electron energy is found from the track momentum

CPP

E = eleCollection[0].getEnergy();

and the second where the electron energy comes from the energy of the associated ECAL cluster

CPP

E = eleCollection[0].getClusters()[0].getEnergy();

The benchmark plots for these two scenarios are overlaid on the same canvas. For the larger Q2 file, the two scenarios perform similarly (for March 2025 files considered here) but when looking at the low Q2 file, some differences become apparent. As one might expect, the Double Angle and JB methods are unaffected by this change, as they do not use the scattered electron energy in their calculation. However, the electron method, and to a lesser extent the Sigma methods see an improvement when using the energy value from the ECAL. The takeaway here is to check which approach gives you a better resolution for your analysis - at low Q2 it’s generally better to rely on the calorimeters, as tracking momentum resolutions are worse at shallow angles.

There are many ways in which the reconstruction of the hadronic final state, which has not been discussed much in this tutorial, could be improved. The reader is invited to look through the current HFS reconstruction code, which constructs the HFS as the sum of particles in the ReconstructedParticles branch, excluding the scattered electron. Note that this code uses boosts to correct for the crossing angle at ePIC - this is something to keep in mind if you are reconstructing the HFS manually, as it strongly impacts the HFS inputs to the reconstruction methods (Pt and E-pz sum).

Challenge

Exercise

Try your own reconstruction of the hadronic final state, for example using a particle-flow algorithm or kinematic fitting of exclusive final states, and compare the resulting resolutions to the default HFS reconstruction. Remember to correct for the 25mrad crossing angle.

Show me the solution

There is no single canonical answer - this is an open-ended, research-level exercise. A reasonable starting point is to reproduce the default HFS (the momentum sum of ReconstructedParticles excluding the scattered electron, boosted to correct for the crossing angle as in the HadronicFinalState code) and then swap in your own particle selection. Feed the resulting pt_had and sigma_h into the calc_jb_method/calc_da_method/calc_sig_method functions and check whether the (reco-true)/true resolutions improve. Improvements from particle flow or exclusive kinematic fitting will hopefully be the subject of a future tutorial.

Key Points

• Some kinematics may favour a calorimeter based electron energy over a tracking based determination - check which is better for your analysis.