1+ // only inside the ROOT prompt / root -x
2+ #ifdef __CLING__
3+ R__LOAD_LIBRARY (podioDict)
4+ R__LOAD_LIBRARY(podioRootIO)
5+ R__LOAD_LIBRARY(libedm4hepDict)
6+ R__LOAD_LIBRARY(libedm4eicDict)
7+ #endif
8+
9+ #include " podio/Frame.h"
10+ #include " podio/ROOTReader.h"
11+ #include " podio/ROOTWriter.h"
12+
13+ #include < edm4hep/MCParticleCollection.h>
14+
15+ #include < fmt/core.h>
16+ #include < fmt/format.h>
17+
18+ #include < iostream>
19+ #include < string>
20+ #include < vector>
21+ #include < cstdlib>
22+ // Global variables for analysis configuration
23+ int events_limit = -1 ; // -1 means process all events
24+ int total_evt_counter = 0 ;
25+
26+ // Function declarations
27+ void print_usage (const char * program_name);
28+ void process_file (const std::string& filename);
29+ void process_event (const podio::Frame& event, int event_number);
30+
31+
32+
33+ int main (int argc, char * argv[]) {
34+ // Simple argument parsing - collect files and look for -n option
35+ std::vector<std::string> input_files;
36+
37+ for (int i = 1 ; i < argc; ++i) {
38+ std::string arg = argv[i];
39+
40+ if (arg == " -n" && i + 1 < argc) {
41+ events_limit = std::atoi (argv[++i]);
42+ fmt::print (" Event limit set to: {}\n " , events_limit);
43+ } else if (arg == " -h" || arg == " --help" ) {
44+ print_usage (argv[0 ]);
45+ return 0 ;
46+ } else if (arg[0 ] != ' -' ) {
47+ input_files.push_back (arg);
48+ } else {
49+ fmt::print (" Unknown option: {}\n " , arg);
50+ print_usage (argv[0 ]);
51+ return 1 ;
52+ }
53+ }
54+
55+ // Check if we have input files
56+ if (input_files.empty ()) {
57+ fmt::print (" Error: No input files provided\n " );
58+ print_usage (argv[0 ]);
59+ return 1 ;
60+ }
61+
62+ // Process each file
63+ fmt::print (" Processing {} file(s)\n " , input_files.size ());
64+ for (const auto & filename : input_files) {
65+ fmt::print (" \n === Processing file: {} ===\n " , filename);
66+ process_file (filename);
67+
68+ // Check if we've reached the event limit
69+ if (events_limit > 0 && total_evt_counter >= events_limit) {
70+ fmt::print (" \n Reached event limit of {}, stopping.\n " , events_limit);
71+ break ;
72+ }
73+ }
74+
75+ fmt::print (" \n Total events processed: {}\n " , total_evt_counter);
76+ return 0 ;
77+ }
78+
79+ void print_usage (const char * program_name) {
80+ fmt::print (" Usage: {} [options] file1.root [file2.root ...]\n " , program_name);
81+ fmt::print (" Options:\n " );
82+ fmt::print (" -n <number> Process only <number> events (default: all)\n " );
83+ fmt::print (" -h Show this help message\n " );
84+ fmt::print (" \n Note: Options and files can be mixed in any order\n " );
85+ fmt::print (" Example: {} file1.root -n 100 file2.root file3.root\n " , program_name);
86+ }
87+
88+ void process_file (const std::string& filename) {
89+ try {
90+ // Open the ROOT file
91+ auto tfile = TFile::Open (filename.c_str ());
92+ tfile->Print ();
93+ fmt::print (" =====here=======\n " );
94+
95+ auto reader = podio::ROOTReader ();
96+ reader.openFile (filename);
97+
98+ const auto nEvents = reader.getEntries (podio::Category::Event);
99+ fmt::print (" File contains {} events\n " , nEvents);
100+
101+ // Process events
102+ for (unsigned i = 0 ; i < nEvents; ++i) {
103+ // Check event limit
104+ if (events_limit > 0 && total_evt_counter >= events_limit) {
105+ break ;
106+ }
107+
108+ auto event = podio::Frame (reader.readNextEntry (podio::Category::Event));
109+
110+ // For the single event we show how to get access to event level parameters
111+ if (i==0 ) {
112+ const auto & cols = event.getParameterKeys <std::string>();
113+ fmt::print (" ===== Collections =====\n " );
114+ for (const auto & col : cols) {
115+ fmt::print (" {} {}\n " , col, event.getParameter <std::string>(col).value_or (" None" ));
116+ }
117+ fmt::print (" =======================\n " );
118+ }
119+ process_event (event, i);
120+ total_evt_counter++;
121+ }
122+ } catch (const std::exception& e) {
123+ fmt::print (" Error processing file {}: {}\n " , filename, e.what ());
124+ }
125+ }
126+
127+ void process_event (const podio::Frame& event, int event_number) {
128+ // Get MCParticles collection
129+ try {
130+ auto & mcParticles = event.get <edm4hep::MCParticleCollection>(" MCParticles" );
131+
132+ if (event_number < 5 ) { // Debug output for first few events
133+ fmt::print (" \n --- Event {} ---\n " , event_number);
134+ fmt::print (" Number of MCParticles: {}\n " , mcParticles.size ());
135+
136+ // Print detailed information for first few particles
137+ int particle_count = 0 ;
138+ for (const auto & particle : mcParticles) {
139+ // if (particle_count >= 5) break; // Only show first 5 particles
140+
141+ fmt::print (" \n Particle {}:\n " , particle_count);
142+ fmt::print (" PDG: {}\n " , particle.getPDG ());
143+ fmt::print (" IDX: {}\n " , particle.getObjectID ().index );
144+ fmt::print (" Charge: {:.3f}\n " , particle.getCharge ());
145+ fmt::print (" Mass: {:.6f} GeV\n " , particle.getMass ());
146+ fmt::print (" GeneratorStatus: {}\n " , particle.getGeneratorStatus ());
147+ fmt::print (" SimulatorStatus: {}\n " , particle.getSimulatorStatus ());
148+
149+ // Momentum
150+ auto momentum = particle.getMomentum ();
151+ fmt::print (" Momentum: ({:.3f}, {:.3f}, {:.3f}) GeV\n " ,
152+ momentum.x , momentum.y , momentum.z );
153+
154+ // Vertex
155+ auto vertex = particle.getVertex ();
156+ fmt::print (" Vertex: ({:.3f}, {:.3f}, {:.3f}) mm\n " ,
157+ vertex.x , vertex.y , vertex.z );
158+
159+ // Endpoint
160+ auto endpoint = particle.getEndpoint ();
161+ fmt::print (" Endpoint: ({:.3f}, {:.3f}, {:.3f}) mm\n " ,
162+ endpoint.x , endpoint.y , endpoint.z );
163+
164+ // Time
165+ fmt::print (" Time: {:.3f} ns\n " , particle.getTime ());
166+
167+ // Parent/daughter relationships
168+ fmt::print (" Number of parents: {}\n " , particle.getParents ().size ());
169+ fmt::print (" Number of daughters: {}\n " , particle.getDaughters ().size ());
170+
171+ if (particle.getDaughters ().size () > 1 ) {
172+ for (const auto & daughter : particle.getDaughters ())
173+ fmt::print (" daughter: index: {} PDG: {}\n " , daughter.getObjectID ().index , daughter.getPDG ());
174+ }
175+
176+ particle_count++;
177+ }
178+ }
179+
180+ // TODO: Add your analysis logic here
181+ // For example, finding Lambda particles, analyzing decays, etc.
182+
183+ } catch (const std::exception& e) {
184+ fmt::print (" Error accessing MCParticles in event {}: {}\n " , event_number, e.what ());
185+ }
186+ }
187+
188+ // ---------------------------------------------------------------------------
189+ // ROOT-macro entry point.
190+ // Call it from the prompt: root -x -l -b -q 'test_edm4eic.cxx("file.root",100)'
191+ // ---------------------------------------------------------------------------
192+ void csv_mcdis (const char * infile, int events = -1 ) {
193+
194+ fmt::print (" 'test_edm4eic' entry point is used\n " );
195+ fmt::print (" infile: {}" , infile);
196+ fmt::print (" events: {}" , events);
197+
198+ events_limit = events; // reuse the global controls
199+ total_evt_counter = 0 ;
200+
201+ process_file (infile);
202+
203+ fmt::print (" \n Total events processed: {}\n " , total_evt_counter);
204+ }
0 commit comments