Skip to content

Commit feeee53

Browse files
DraTeotsDraTeots
committed
EDM4EIC CPP analysis
1 parent 8a4a6b7 commit feeee53

File tree

4 files changed

+242
-0
lines changed

4 files changed

+242
-0
lines changed

‎analysis/cpp-example/e01_edm4hep.cxx‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@
44
R__LOAD_LIBRARY(podioDict)
55
R__LOAD_LIBRARY(podioRootIO)
66
R__LOAD_LIBRARY(libedm4hepDict)
7+
R__LOAD_LIBRARY(libedm4eicDict)
78
#include "podio/EventStore.h"
89
#include "podio/ROOTReader.h"
910
#include "edm4hep/utils/kinematics.h"

‎analysis/cpp-example/test_4eic.cxx‎

Lines changed: 147 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,147 @@
1+
// singularity exec -B /volatile/eic/romanov/meson-structure-2025-06:/volatile/eic/romanov/meson-structure-2025-06 /cvmfs/singularity.opensciencegrid.org/eicweb/eic_xl:nightly bash
2+
// source /opt/detector/epic-main/bin/thisepic.sh
3+
// root -x -l -b -q e01_edm4hep.cxx\(\"/volatile/eic/romanov/meson-structure-2025-06/analysis/cpp-example/data/DRICHPID.root\"\)
4+
R__LOAD_LIBRARY(podioDict)
5+
R__LOAD_LIBRARY(podioRootIO)
6+
R__LOAD_LIBRARY(libedm4hepDict)
7+
R__LOAD_LIBRARY(libedm4eicDict)
8+
#include "podio/EventStore.h"
9+
#include "podio/ROOTReader.h"
10+
#include "edm4hep/utils/kinematics.h"
11+
12+
// #define _AEROGEL_
13+
14+
#define _NPE_REFERENCE_ 211
15+
//#define _NPE_REFERENCE_ (11)
16+
//#define _NPE_REFERENCE_ 321
17+
18+
void e01_edm4hep(const char *ifname, const char *ofname = nullptr)
19+
{
20+
// open recon output .root file with podio
21+
podio::EventStore podio_store;
22+
podio::ROOTReader podio_reader;
23+
podio_reader.openFile(ifname);
24+
podio_store.setReader(&podio_reader);
25+
26+
auto np = new TH1D("np", "Photon count", 50, 0, 50);
27+
28+
29+
// event loop
30+
for(unsigned ev=0; ev<podio_reader.getEntries(); ev++) {
31+
if(ev%100==0) printf("read event %d\n",ev);
32+
33+
// get collections
34+
auto& cherenkovs = podio_store.get<edm4eic::CherenkovParticleIDCollection>("DRICHPID");
35+
auto& mctracks = podio_store.get<edm4hep::MCParticleCollection>("MCParticles");
36+
37+
// Then the Cherenkov-to-reconstructed mapping;
38+
// FIXME: may want to use Cherenkov-to-simulated mapping to start with, for the debugging purposes;
39+
// FIXME: if we loop over reconstructed tracks, rather than MC particles, then we
40+
// have the 1-1 relation edm4hep::ReconstructedParticle::getParticleIDUsed(),
41+
// which returns the type edm4hep::ParticleID
42+
std::map<edm4hep::MCParticle,edm4eic::CherenkovParticleID> rc2cherenkov;
43+
for(const auto &pid : cherenkovs)
44+
rc2cherenkov[pid.getAssociatedParticle()] = pid;
45+
46+
// Loop through all MC tracks;
47+
for(const auto &mctrack : mctracks) {
48+
// FIXME: consider only primaries for now?; equivalent to mctrack.getGeneratorStatus()==1?
49+
if (mctrack.parents_size()>0) continue;
50+
51+
edm4eic::CherenkovParticleID cherenkov;
52+
if(rc2cherenkov.find(mctrack) != rc2cherenkov.end()) cherenkov = rc2cherenkov[mctrack];
53+
else continue;
54+
55+
double pp = edm4hep::utils::p(mctrack);
56+
double m = mctrack.getMass();
57+
58+
//printf("m=%5.3f p=%5.1f (%4d) \n", m, pp, mctrack.getPDG());
59+
60+
// Loop through all of the mass hypotheses available for this reconstructed track;
61+
{
62+
const edm4eic::CherenkovPdgHypothesis *best = 0;
63+
64+
for(const auto &option : cherenkov.getOptions()) {
65+
66+
if (option.radiator != id) continue;
67+
68+
// Skip electron hypothesis; of no interest here;
69+
//if (abs(option.pdg) == 11) continue;
70+
71+
if (abs(option.pdg) == _NPE_REFERENCE_) {
72+
np->Fill(option.npe);
73+
74+
if (ofname) npvector.push_back(option.npe);
75+
} //if
76+
77+
if (!best || option.weight > best->weight) best = &option;
78+
printf("radiator %3d (pdg %5d): weight %7.2f, npe %7.2f\n",
79+
option.radiator, option.pdg, option.weight, option.npe);
80+
} //for option
81+
printf("\n");
82+
83+
// Check whether the true PDG got a highest score;
84+
if (!best || best->pdg != mctrack.getPDG()) false_assignment_stat[best->npe >= 5 ? 0 : 1]++;
85+
86+
// This assumes of course that at least one radiator was requested in juggler;
87+
double rindex = cherenkov.getAngles()[id].rindex;
88+
double theta = cherenkov.getAngles()[id].theta;
89+
double lambda = cherenkov.getAngles()[id].wavelength;
90+
double argument = sqrt(pp*pp + m*m)/(rindex*pp);
91+
double thp = fabs(argument) <= 1.0 ? acos(argument) : theta;
92+
93+
th->Fill(1000 * theta);
94+
//if (mctrack.getPDG() == 321)
95+
dt->Fill(1000* (theta - thp));
96+
ri->Fill(rindex - 1.0);
97+
wl->Fill(lambda);//rindex - 1.0);
98+
printf("<n> ~ %8.6f, <th> = %7.2f [mrad]\n", rindex - 1.0, 1000*thp);
99+
100+
if (ofname) thvector.push_back(theta - thp);
101+
}
102+
} //for track
103+
104+
// next event
105+
podio_store.clear();
106+
podio_reader.endOfEvent();
107+
} //for ev
108+
109+
// end of event loop
110+
printf("%3d (%3d) false out of %d\n", false_assignment_stat[0], false_assignment_stat[1], podio_reader.getEntries());
111+
podio_reader.closeFile();
112+
113+
// write
114+
if (ofname) {
115+
116+
auto *ofdata = new TFile(ofname, "RECREATE");
117+
118+
if (!ofdata) {
119+
printf("was not able to create output file '%s'\n", ofname);
120+
exit(0);
121+
} //if
122+
auto *ot = new TTree("t", "My tree");
123+
124+
double thbff, npbff;
125+
ot->Branch("th", &thbff, "th/D");
126+
ot->Branch("np", &npbff, "np/D");
127+
128+
for(unsigned iq=0; iq<thvector.size(); iq++) {
129+
thbff = thvector[iq];
130+
npbff = npvector[iq];
131+
132+
ot->Fill();
133+
} //for iq
134+
135+
ot->Write();
136+
ofdata->Close();
137+
exit(0);
138+
} else {
139+
auto cv = new TCanvas("cv", "", 1700, 500);
140+
cv->Divide(5, 1);
141+
cv->cd(1); np->Draw(); np->Fit("gaus");
142+
cv->cd(2); th->Draw("SAME"); th->Fit("gaus");
143+
cv->cd(3); ri->Draw(); ri->Fit("gaus");
144+
cv->cd(4); dt->Draw(); dt->Fit("gaus");
145+
cv->cd(5); wl->Draw(); //dt->Fit("gaus");
146+
} //if
147+
} // evaluation()

‎analysis/zdc-lambda/CMakeLists.txt‎

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,29 @@
1+
# ------------- CMakeLists.txt -----------------
2+
cmake_minimum_required(VERSION 3.18)
3+
project(cpp_meson_structure LANGUAGES CXX)
4+
5+
# ----- compiler settings ------------------------------------------------------
6+
set(CMAKE_CXX_STANDARD 20)
7+
set(CMAKE_CXX_STANDARD_REQUIRED ON)
8+
set(CMAKE_CXX_EXTENSIONS OFF)
9+
10+
# -----------------------------------------------------------------------------#
11+
# Dependencies #
12+
# -----------------------------------------------------------------------------#
13+
14+
find_package(ROOT REQUIRED COMPONENTS Core RIO Tree Hist Gpad Graf)
15+
find_package(podio CONFIG REQUIRED)
16+
find_package(EDM4HEP CONFIG REQUIRED)
17+
find_package(EDM4EIC CONFIG REQUIRED)
18+
19+
# -----------------------------------------------------------------------------#
20+
# Build executable #
21+
# -----------------------------------------------------------------------------#
22+
add_executable(e01_edm4hep test_edm4eic.cxx)
23+
24+
target_link_libraries(e01_edm4hep
25+
PRIVATE
26+
ROOT::Core ROOT::RIO ROOT::Tree ROOT::Hist ROOT::Gpad
27+
podio::podio podio::podioRootIO
28+
EDM4HEP::edm4hep
29+
EDM4EIC::edm4eic)

‎analysis/zdc-lambda/test_edm4eic.cxx‎

Lines changed: 65 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,65 @@
1+
#include "podio/Frame.h"
2+
#include "podio/ROOTReader.h"
3+
#include "podio/ROOTWriter.h"
4+
5+
#include <edm4hep/MCParticleCollection.h>
6+
7+
#include <iostream>
8+
#include <string>
9+
10+
void readCollection() {
11+
// Start reading the input
12+
auto reader = podio::ROOTReader();
13+
reader.openFile("associations.root");
14+
15+
const auto nEvents = reader.getEntries(podio::Category::Event);
16+
17+
for (unsigned i = 0; i < 5; ++i) {
18+
auto store = podio::Frame(reader.readNextEntry(podio::Category::Event));
19+
20+
auto& mcParticles = store.get<edm4hep::MCParticleCollection>("MCParticles");
21+
for (const auto& particle : mcParticles) {
22+
std::cout << particle.getPDG() << std::endl;
23+
}
24+
// if (clusters.isValid()) {
25+
// for (const auto& cluster : clusters) {
26+
// if (cluster.isAvailable()) {
27+
// for (const auto& hit : cluster.Hits()) {
28+
// if (hit.isAvailable()) {
29+
// throw std::runtime_error("Hit is available, although it has not been written");
30+
// }
31+
// }
32+
// }
33+
// }
34+
// } else {
35+
// throw std::runtime_error("Collection 'clusters' should be present");
36+
// }
37+
//
38+
// // Test for subset collections
39+
// auto& hits_subset = store.get<ExampleHitCollection>("hits_subset");
40+
// if (hits_subset.isValid()) {
41+
// if (!hits_subset.isSubsetCollection()) {
42+
// throw std::runtime_error("hits_subset should be a subset collection");
43+
// }
44+
//
45+
// if (hits_subset.size() != 2) {
46+
// throw std::runtime_error("subset collection should have original size");
47+
// }
48+
//
49+
// for (const auto& hit : hits_subset) {
50+
// if (hit.isAvailable()) {
51+
// throw std::runtime_error("Hit is available, although it has not been written");
52+
// }
53+
// }
54+
// } else {
55+
// throw std::runtime_error("Collection 'hits_subset' should be present");
56+
// }
57+
}
58+
}
59+
60+
int main() {
61+
62+
readCollection();
63+
64+
return 0;
65+
}

0 commit comments

Comments
 (0)