Skip to content
Snippets Groups Projects

Part5 analysis docs

Merged Part5 analysis docs
wdconinc-master-patch-73409 into master
Merged Wouter Deconinck requested to merge wdconinc-master-patch-73409 into master
Compare
and
1 file
+ 304
− 0
Compare changes
  • Side-by-side
  • Inline
src/docs/part5/reconstruction_analysis.md 0 → 100644
+ 304
− 0
---
title: 'Reconstruction Analysis'
---
This part of the tutorial outlines how reconstruction output from full simulations can be used for analysis by physics working groups.
At the end of this tutorial you should able to:
- determine the location of relevant full simulation reconstruction files,
- understand the file structure of the full simulation reconstruction output,
- use one of several approaches to perform low-level analysis on the reconstruction output.
## Location of full simulation reconstruction output
Full simulation outputs are stored on one of several locations:
- S3: https://dtn01.sdcc.bnl.gov:9000/minio/eictest/ATHENA/ (web interface)
- xrootd: root://sci-xrootd.jlab.org//osgpool/eic/ATHENA/ (no preview currently available)
Data at these location is organized in a predicatable structure:
- `EVGEN/` contains all initial generated events (in hepmc3 or other format),
- `FULL/` contains the raw full simulation output without any reconstruction,
- `RECO/` contains the reconstruction output.
Under the different top-level directories, you can find the different physics processes:
- `SINGLE/` contains single particle initial states,
- `DIS/` contains DIS events.
For details on accessing these locations, please refer to https://doc.athena-eic.org/en/latest/howto/.
For the purpose of this tutorial we will use the S3 interface since it does not rely on a mirroring process at Jefferson Lab.
## File structure of the full simulation reconstruction output
Each reconstruction output file has essentially the same structure, defined by the EIC Data Model. This structure can be retrieved with `rootls`, e.g.
```console
root -l s3https://dtn01.sdcc.bnl.gov:9000/eictest/ATHENA/RECO/SINGLE/pi+/1GeV/45to135deg/pi+_1GeV_45to135deg.0001.root
```
which produces the following output, two trees:
```console
events metadata
```
The `events` tree is of course what we are interested in. We can explore its top-level structure as follows. We first start a ROOT session:
```console
root -l s3https://dtn01.sdcc.bnl.gov:9000/eictest/ATHENA/RECO/SINGLE/pi+/1GeV/45to135deg/pi+_1GeV_45to135deg.0001.root
```
and then list the top-level branches in the `events` tree:
```console
events->Print("toponly")
```
This will display a large number of branches (and their size):
```console
******************************************************************************
*Tree :events : Events tree *
*Entries : 1000002 : Total = 34023140832 bytes File Size = 8907188433 *
* : : Tree compression factor = 3.82 *
******************************************************************************
branch: mcparticles2 471802624
branch: TrackerBarrelHits2 184032064
branch: ReconstructedParticles 2106267
branch: ReconstructedParticles#0 1447954
branch: ReconstructedParticles#1 1447954
branch: ReconstructedParticles#2 1447954
branch: ReconstructedParticles#3 1447954
branch: EcalBarrelHitsSimpleDigi 93024680
branch: EcalBarrelHitsSimpleReco 142564938
branch: EcalBarrelHitsSimpleReco#0 6351262
branch: EcalEndcapNHitsDigi 17873069
branch: EcalEndcapNHitsReco 9750180
branch: EcalEndcapNHitsReco#0 2294175
branch: EcalEndcapNClusterHits 3686937
branch: EcalEndcapNClusterHits#0 1594591
branch: EcalEndcapNClusters 4128541
branch: EcalEndcapNClusters#0 1444440
branch: EcalEndcapNClusters#1 1444440
branch: EcalEndcapPHitsDigi 43298624
branch: EcalEndcapPHitsReco 51470328
branch: EcalEndcapPHitsReco#0 3366251
branch: EcalEndcapPHitsRecoXY 48939215
branch: EcalEndcapPHitsRecoXY#0 3383190
branch: EcalEndcapPClusterHits 6539903
branch: EcalEndcapPClusterHits#0 1633199
branch: EcalEndcapPClusters 4494705
branch: EcalEndcapPClusters#0 1444440
branch: EcalEndcapPClusters#1 1444440
branch: EcalBarrelHitsDigi 68681726
branch: EcalBarrelHitsReco 441682325
branch: EcalBarrelClusterHits 22471306
branch: EcalBarrelLayers 7351806
branch: EcalBarrelLayers#0 1441221
branch: EcalBarrelClusters 7131448
branch: EcalBarrelClusters#0 1443270
branch: EcalBarrelClusters#1 1443270
branch: EcalBarrelScFiHitsDigi 1636334484
branch: EcalBarrelScFiHitsReco 3932534929
branch: EcalBarrelScFiHitsReco#0 14507111
branch: EcalBarrelScFiGridReco 815861620
branch: EcalBarrelScFiGridReco#0 9424066
branch: EcalBarrelScFiClusterHits 396488274
branch: EcalBarrelScFiClusterHits#0 6066163
branch: EcalBarrelScFiClusters 41044730
branch: EcalBarrelScFiClusters#0 1447954
branch: EcalBarrelScFiClusters#1 1447954
branch: HcalBarrelHitsDigi 13328585
branch: HcalBarrelHitsReco 15260047
branch: HcalBarrelHitsReco#0 2425669
branch: HcalBarrelHitsRecoXY 13434597
branch: HcalBarrelHitsRecoXY#0 2410125
branch: HcalBarrelClusterHits 2352864
branch: HcalBarrelClusterHits#0 1451387
branch: HcalBarrelClusters 2687875
branch: HcalBarrelClusters#0 1443270
branch: HcalBarrelClusters#1 1443270
branch: HcalElectronEndcapHitsDigi 5707858
branch: HcalElectronEndcapHitsReco 8514858
branch: HcalElectronEndcapHitsReco#0 1928251
branch: HcalElectronEndcapHitsRecoXY 7773148
branch: HcalElectronEndcapHitsRecoXY#0 1922562
branch: HcalElectronEndcapClusterHits 2375296
branch: HcalElectronEndcapClusterHits#0 1459130
branch: HcalElectronEndcapClusters 2726638
branch: HcalElectronEndcapClusters#0 1452637
branch: HcalElectronEndcapClusters#1 1452637
branch: HcalHadronEndcapHitsDigi 4112230
branch: HcalHadronEndcapHitsReco 5003312
branch: HcalHadronEndcapHitsReco#0 1689929
branch: HcalHadronEndcapHitsRecoXY 4723184
branch: HcalHadronEndcapHitsRecoXY#0 1689633
branch: HcalHadronEndcapClusterHits 2328711
branch: HcalHadronEndcapClusterHits#0 1454456
branch: HcalHadronEndcapClusters 2690491
branch: HcalHadronEndcapClusters#0 1450880
branch: HcalHadronEndcapClusters#1 1450880
branch: TrackerBarrelRawHits 24139635
branch: TrackerEndcapRawHits 3166129
branch: VertexBarrelRawHits 29242781
branch: VertexEndcapRawHits 2219822
branch: TrackerBarrelRecHits 64929921
branch: TrackerEndcapRecHits 6218655
branch: VertexBarrelRecHits 75481339
branch: VertexEndcapRecHits 3798116
branch: ReconstructedParticlesInitFromTruth 1958375
branch: outputTrackParameters 2035687
branch: ForwardRICHHits2 20341908
branch: ForwardRICHHitsDigi 1500666
branch: ForwardRICHHitsReco 1858812
branch: EcalEndcapNHits 2255311
branch: EcalEndcapPHits 2255311
branch: EcalBarrelHits 2248873
branch: EcalBarrelScFiHits 2272009
branch: HcalBarrelHits 2248873
branch: HcalHadronEndcapHits 2280802
branch: HcalElectronEndcapHits 2291933
branch: TrackerEndcapHits 2558412
branch: TrackerBarrelHits 2558412
branch: VertexBarrelHits 2554303
branch: VertexEndcapHits 2554303
branch: ForwardRICHHits 2544344
```
During the development of the reconstruction, there are more branches enabled here than are strictly necessary (e.g. simulated and digitized hits, intermediate reconstruction parameters). These are all available for analysis (with fixed interfaces). In this tutorial we will focus on a few branches in particular:
- ReconstructedParticles: contains the results from track finding and fitting,
- EcalBarrelClusters: contains the results from the barrel Ecal cluster finding,
- EcalBarrelScFiClusters: contains the results from the barrel Ecal ScFi cluster finding.
We can inspect each of these three branches in more detail (some information removed for formatting)
```console
root [14] events->Print("ReconstructedParticles*")
******************************************************************************
*Tree :events : Events tree *
*Entries : 1000002 : Total = 34023140832 bytes File Size = 8907188433 *
* : : Tree compression factor = 3.82 *
******************************************************************************
*Br 0 :ReconstructedParticles : Int_t ReconstructedParticles_ *
*Br 1 :ReconstructedParticles.pid : Long64_t pid[ReconstructedParticles_] *
*Br 2 :ReconstructedParticles.energy : *
* | Double_t energy[ReconstructedParticles_] *
*Br 3 :ReconstructedParticles.p.x : Double_t x[ReconstructedParticles_] *
*Br 4 :ReconstructedParticles.p.y : Double_t y[ReconstructedParticles_] *
*Br 5 :ReconstructedParticles.p.z : Double_t z[ReconstructedParticles_] *
*Br 6 :ReconstructedParticles.charge : *
* | Double_t charge[ReconstructedParticles_] *
*Br 7 :ReconstructedParticles.mass : *
* | Double_t mass[ReconstructedParticles_] *
*Br 8 :ReconstructedParticles.clusters_begin : *
* | UInt_t clusters_begin[ReconstructedParticles_] *
*Br 9 :ReconstructedParticles.clusters_end : *
* | UInt_t clusters_end[ReconstructedParticles_] *
*Br 10 :ReconstructedParticles.tracks_begin : *
* | UInt_t tracks_begin[ReconstructedParticles_] *
*Br 11 :ReconstructedParticles.tracks_end : *
* | UInt_t tracks_end[ReconstructedParticles_] *
*Br 12 :ReconstructedParticles.particles_begin : *
* | UInt_t particles_begin[ReconstructedParticles_] *
*Br 13 :ReconstructedParticles.particles_end : *
* | UInt_t particles_end[ReconstructedParticles_] *
*............................................................................*
root [6] events->Print("EcalBarrelClusters*")
******************************************************************************
*Tree :events : Events tree *
*Entries : 500002 : Total = 22026619265 bytes File Size = 6339571464 *
* : : Tree compression factor = 3.47 *
******************************************************************************
*Br 0 :EcalBarrelClusters : Int_t EcalBarrelClusters_ *
*Br 1 :EcalBarrelClusters.clusterID : Int_t clusterID[EcalBarrelClusters_]*
*Br 2 :EcalBarrelClusters.nhits : Int_t nhits[EcalBarrelClusters_] *
*Br 3 :EcalBarrelClusters.energy : Double_t energy[EcalBarrelClusters_] *
*Br 4 :EcalBarrelClusters.edep : Double_t edep[EcalBarrelClusters_] *
*Br 5 :EcalBarrelClusters.radius : Double_t radius[EcalBarrelClusters_] *
*Br 6 :EcalBarrelClusters.skewness : *
* | Double_t skewness[EcalBarrelClusters_] *
*Br 7 :EcalBarrelClusters.leakcorr : *
* | Double_t leakcorr[EcalBarrelClusters_] *
*Br 8 :EcalBarrelClusters.eta : Double_t eta[EcalBarrelClusters_] *
*Br 9 :EcalBarrelClusters.position.x : Double_t x[EcalBarrelClusters_] *
*Br 10 :EcalBarrelClusters.position.y : Double_t y[EcalBarrelClusters_] *
*Br 11 :EcalBarrelClusters.position.z : Double_t z[EcalBarrelClusters_] *
*Br 12 :EcalBarrelClusters.polar.r : Double_t r[EcalBarrelClusters_] *
*Br 13 :EcalBarrelClusters.polar.theta : *
* | Double_t theta[EcalBarrelClusters_] *
*Br 14 :EcalBarrelClusters.polar.phi : Double_t phi[EcalBarrelClusters_] *
*Br 15 :EcalBarrelClusters.cl_theta : *
* | Double_t cl_theta[EcalBarrelClusters_] *
*Br 16 :EcalBarrelClusters.cl_phi : Double_t cl_phi[EcalBarrelClusters_] *
*Br 17 :EcalBarrelClusters.hits_begin : *
* | UInt_t hits_begin[EcalBarrelClusters_] *
*Br 18 :EcalBarrelClusters.hits_end : UInt_t hits_end[EcalBarrelClusters_] *
*Br 19 :EcalBarrelClusters.layers_begin : *
* | UInt_t layers_begin[EcalBarrelClusters_] *
*Br 20 :EcalBarrelClusters.layers_end : *
* | UInt_t layers_end[EcalBarrelClusters_] *
*............................................................................*
```
The ReconstructedParticles branch contains the momentum, e.g. `ReconstructedParticles.p.x` for the x-component, and references to other entities related to this particle (clusters, tracks, particles).
Note: Due to conditions that are currently being addressed, some momentum vectors use `p.x` whereas others use `p.px`.
## Analysis of full simulation reconstruction output with traditional ROOT commands
After opening the reconstruction output file, let's make some pretty plots. We start with:
```console
root -l s3https://dtn01.sdcc.bnl.gov:9000/eictest/ATHENA/RECO/SINGLE/pi+/1GeV/45to135deg/pi+_1GeV_45to135deg.0001.root
```
and in the ROOT session:
```console
root [1] .ls
TNetXNGFile** root://sci-xrootd.jlab.org//osgpool/eic/ATHENA/RECO/SINGLE/pi+/1GeV/45to135deg/pi+_1GeV_45to135deg.0001.root data file
TNetXNGFile* root://sci-xrootd.jlab.org//osgpool/eic/ATHENA/RECO/SINGLE/pi+/1GeV/45to135deg/pi+_1GeV_45to135deg.0001.root data file
KEY: TTree events;1 Events tree
KEY: TTree metadata;1 Metadata tree
events->Draw("ReconstructedParticlesInitFromTruth.p.px")
events->Draw("EcalBarrelClusters.polar.phi:EcalBarrelClusters.polar.theta", "EcalBarrelClusters.edep", "colz")
```
These types of plots will likely be limited to simple data inspection.
## Analysis of full simulation reconstruction output with RDataFrame commands
For a more advanced analysis, you can take advantage of the RDataFrame features, such as in this (shortened) DIS example. The [original](https://eicweb.phy.anl.gov/EIC/benchmarks/physics_benchmarks/-/blob/master/benchmarks/dis/analysis/dis_electrons.cxx) is used in our CI system and determines DIS parameters for every change to the detector geometry, simulation, digitization, or reconstruction.
```console
auto momenta_from_reconstruction(const std::vector<eic::ReconstructedParticleData>& parts) {
std::vector<ROOT::Math::PxPyPzEVector> momenta{parts.size()};
std::transform(parts.begin(), parts.end(), momenta.begin(), [](const auto& part) {
return ROOT::Math::PxPyPzEVector{part.p.x, part.p.y, part.p.z, part.energy};
});
return momenta;
}
auto convertMtoE(const std::vector<ROOT::Math::PxPyPzMVector>& mom) {
std::vector<ROOT::Math::PxPyPzEVector> momenta{mom.size()};
std::transform(mom.begin(), mom.end(), momenta.begin(), [](const auto& part) {
return ROOT::Math::PxPyPzEVector{part.Px(), part.Py(), part.Pz(), part.energy()};
});
return momenta;
}
bool sort_mom_bool(ROOT::Math::PxPyPzEVector &mom1, ROOT::Math::PxPyPzEVector &mom2) {
return mom1.energy() > mom2.energy();
}
auto sort_momenta(const std::vector<ROOT::Math::PxPyPzEVector>& mom) {
std::vector <ROOT::Math::PxPyPzEVector> sort_mom = mom;
sort(sort_mom.begin(), sort_mom.end(), sort_mom_bool);
return sort_mom;
}
auto Q2(const std::vector<ROOT::Math::PxPyPzEVector>& mom) {
std::vector<double> Q2Vec(mom.size() );
ROOT::Math::PxPyPzEVector beamMom = {0, 0, 5, 5};
std::transform(mom.begin(), mom.end(), Q2Vec.begin(), [beamMom](const auto& part) {
return -(part - beamMom).M2();
});
return Q2Vec;
}
ROOT::RDataFrame d("events", "s3https://dtn01.sdcc.bnl.gov:9000/eictest/ATHENA/RECO//DIS/crossDivNrgCrab/DIS_NC_Q2gt10_crossDivNrgCrab_25mRad_5x41_v2.root");
auto d0 = d.Define("p", momenta_from_reconstruction, {"ReconstructedParticles"}).Define("Q2", Q2, {"p"});
auto h_Q2_sim = d0.Histo1D({"h_Q2_sim", "; GeV; counts", 100, -5, 25}, "Q2");
auto& h1_Q2_sim = *h_Q2_sim;
h1_Q2_sim.DrawClone("hist");
```