This lesson is being piloted (Beta version)

EIC Tutorial: Kinematic Reconstruction

Introduction

Overview

Teaching: 15 min
Exercises: 5 min
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:

./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:

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

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 March 2025 campaign: a low Q2 and a high Q2 Neutral Current DIS file.

The software used for EIC simulation/reconstruction/analysis is ever-changing, and an April 2025 update has resulted in some incompatibility between recent builds of eic-shell and simulation files produced before April 2025. We will be working with files from the March 2025 campaign, so to avoid this incompatibility we can start up an older version of eic-shell (after exiting our current eic-shell)

./eic-shell -v 25.03.0-stable

and we can grab our files using -

xrdcp root://dtn-eic.jlab.org//volatile/eic/EPIC/RECO/25.03.1/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/25.03.1/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.

If you download a more recent file, the current build of eic-shell should work.

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:

root -l pythia8NCDIS_18x275_minQ2=1_beamEffects_xAngle=-0.025_hiDiv_1.0001.eicrecon.tree.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:

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

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.

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

Performance Benchmarks

Overview

Teaching: 15 min
Exercises: 5 min
Questions

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

Objectives

  • Produce plots that benchmark the performance of different reconstruction methods

Using ROOTFrameReader 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 here) as an example:

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

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

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

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 ROOTFrameReader would look somthing like this

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

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

for (size_t i = 0; i < reader.getEntries("events"); i++) {
    const auto event = podio::Frame(reader.readNextEntry("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

// PODIO
#include "podio/Frame.h"
#include "podio/ROOTFrameReader.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) {

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

  auto reader = podio::ROOTFrameReader();
  reader.openFiles(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);

  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 = podio::Frame(reader.readNextEntry("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

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

or as

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

// 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.

Key Points

  • Use ROOTFrameReader to process simulation files using the data types implemented in edm4hep/edm4eic.

Manual Reconstruction

Overview

Teaching: 15 min
Exercises: 5 min
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

// PODIO
#include "podio/Frame.h"
#include "podio/ROOTFrameReader.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::ROOTFrameReader();
  reader.openFiles(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 = podio::Frame(reader.readNextEntry("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

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

// 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 branch.

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.

A challenge for the reader: try implementing your own version of the electron method using four vectors - verify that it matches the output of the InclusiveKinematicsElectron branch.

Key Points

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

Optimise Reconstruction

Overview

Teaching: 15 min
Exercises: 5 min
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

// PODIO
#include "podio/Frame.h"
#include "podio/ROOTFrameReader.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::ROOTFrameReader();
  reader.openFiles(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 = podio::Frame(reader.readNextEntry("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");

    auto& ecalClusters = event.get<edm4eic::ClusterCollection>("EcalClusters");

    // 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

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

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

E = eleCollection[0].getEnergy();

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

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). Regardless of the difficulties, it is highly recommended to try your own HFS reconstruction, as there are many improvements to be offered by e.g. particle flow algorithms, or kinematic fitting of exclusive final states, that 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