From 15093eceb2499bf49651b25b9ad13748f92cf1bb Mon Sep 17 00:00:00 2001 From: Matthew Gignac Date: Tue, 7 Apr 2026 15:21:38 -0700 Subject: [PATCH] Revert "Merge pull request #446 from sarahgaiser/add_mother_update" This reverts commit 3763e24306009733145c5f8708e78cc5ddc97c37, reversing changes made to 26a8d8785dc81ecaa80f83713735f7af02f60e35. --- .../stdhep-tools/src/add_mother_full_truth.cc | 659 ++++-------------- 1 file changed, 151 insertions(+), 508 deletions(-) diff --git a/tools/stdhep-tools/src/add_mother_full_truth.cc b/tools/stdhep-tools/src/add_mother_full_truth.cc index cc5aa44c2..258ed6760 100644 --- a/tools/stdhep-tools/src/add_mother_full_truth.cc +++ b/tools/stdhep-tools/src/add_mother_full_truth.cc @@ -6,529 +6,172 @@ #include #include -#include using namespace std; /** - * Adds mother particles (623 scattered electron/muon and 622 A-prime/virtual photon) to events - * and builds proper genealogy from LHE and STDHEP input files. - * Handles both A-prime (nup=7) and radiative (nup=6) cases. + * The tool works for ap and rad MC for fully truth. + * Takes input stdhep file and lhe file, adds new particles and makes them mothers of parentless particles, + * and writes to a new stdhep file. */ +int main(int argc,char** argv) +{ + int nevhep; //!< The event number + vector new_event; -const double ELECTRON_MASS = 0.000511; -const double BEAM_ENERGY = 3.74; + int id_beam = 623; + int id_pair = 622; -struct LHEParticle { - int idhep; - int isthep; - int jmohep[2]; - double phep[5]; // px, py, pz, E, mass -}; + double mass = 0.1; + double energy = 0.1; -enum EventType { - EVENT_APRIME, // nup = 7 - EVENT_RADIATIVE, // nup = 6 - EVENT_UNKNOWN -}; + int c; -bool read_lhe_event_aprime(FILE* lhe_file, int nup, LHEParticle& scattered_electron, - LHEParticle& aprime, int id_pair) { - char line[1000]; - - bool found_aprime = false; - bool found_scattered = false; - - // Read all particles in the event - for (int i = 0; i < nup; i++) { - if (fgets(line, 1000, lhe_file) == NULL) return false; - - LHEParticle temp; - sscanf(line, "%d %d %d %d %*d %*d %lf %lf %lf %lf %lf", - &temp.idhep, &temp.isthep, - &temp.jmohep[0], &temp.jmohep[1], - &temp.phep[0], &temp.phep[1], &temp.phep[2], - &temp.phep[3], &temp.phep[4]); - - // Find A-prime (id_pair, status 2) - if (temp.idhep == id_pair && temp.isthep == 2) { - aprime = temp; - found_aprime = true; - } - - // Find scattered electron (11, status 1, energy < beam energy) - if (temp.idhep == 11 && temp.isthep == 1 && temp.phep[3] < BEAM_ENERGY) { - scattered_electron = temp; - found_scattered = true; - } - } - - if (!found_aprime || !found_scattered) { - fprintf(stderr, "Error: Could not find required particles in A-prime LHE event\n"); - fprintf(stderr, " Found pair particle (%d): %d, Found scattered electron: %d\n", - id_pair, found_aprime, found_scattered); - return false; - } - - return true; -} + while ((c = getopt(argc,argv,"hi1:i2:")) != -1) + switch (c) + { + case 'h': + printf("-h: print this help\n"); + printf("-i: PDG ID of mother\n"); + return(0); + break; + case 'i1': + id_beam = atoi(optarg); + break; + case 'i2': + id_pair = atoi(optarg); + break; + case '?': + printf("Invalid option or missing option argument; -h to list options\n"); + return(1); + default: + abort(); + } -bool read_lhe_event_radiative(FILE* lhe_file, int nup, LHEParticle& scattered_lepton, - LHEParticle& electron, LHEParticle& positron) { - char line[1000]; - - bool found_scattered = false; - bool found_electron = false; - bool found_positron = false; - - // Read all particles in the event - for (int i = 0; i < nup; i++) { - if (fgets(line, 1000, lhe_file) == NULL) return false; - - LHEParticle temp; - sscanf(line, "%d %d %d %d %*d %*d %lf %lf %lf %lf %lf", - &temp.idhep, &temp.isthep, - &temp.jmohep[0], &temp.jmohep[1], - &temp.phep[0], &temp.phep[1], &temp.phep[2], - &temp.phep[3], &temp.phep[4]); - - // Find scattered lepton (13, status 1, energy < beam energy) - if (temp.idhep == 13 && temp.isthep == 1 && temp.phep[3] < BEAM_ENERGY) { - scattered_lepton = temp; - found_scattered = true; - } - - // Find electron (11, status 1) - if (temp.idhep == 11 && temp.isthep == 1) { - electron = temp; - found_electron = true; - } - - // Find positron (-11, status 1) - if (temp.idhep == -11 && temp.isthep == 1) { - positron = temp; - found_positron = true; - } - } - - if (!found_scattered || !found_electron || !found_positron) { - fprintf(stderr, "Error: Could not find required particles in radiative LHE event\n"); - fprintf(stderr, " Found scattered lepton: %d, Found electron: %d, Found positron: %d\n", - found_scattered, found_electron, found_positron); - return false; - } - - return true; -} + if (argc-optind < 3) + { + printf(" \n"); + return 1; + } + + int n_events; + int istream = 0; + int ostream = 1; + + n_events = open_read(argv[optind], istream); + + open_write(argv[optind+2], ostream, n_events); + + FILE * in_file = fopen(argv[optind+1], "r"); + + while (true) { + //! read event from stdhep file + if (!read_next(istream)) { + close_read(istream); + fclose(in_file); + close_write(ostream); + return(0); + } -bool read_lhe_event(FILE* lhe_file, EventType& event_type, int& nup, - LHEParticle& scattered_particle, LHEParticle& pair_or_electron, - LHEParticle& positron, int id_pair) { + /** + * build event + * entry 1: mother particle 1 + * entry 2: mother particle 2 + * other entries: copy of input stdhep event + */ + struct stdhep_entry *temp1 = new struct stdhep_entry; + temp1->isthep = 3; //!> documentation particle + temp1->idhep = id_beam; + for (int j = 0; j < 2; j++) temp1->jmohep[j] = 0; + for (int j = 0; j < 2; j++) temp1->jdahep[j] = 0; + for (int j = 0; j < 5; j++) temp1->phep[j] = 0.0; + temp1->phep[3] += energy; + temp1->phep[4] += mass; + for (int j = 0; j < 4; j++) temp1->vhep[j] = 0.0; + new_event.push_back(*temp1); + + struct stdhep_entry *temp2 = new struct stdhep_entry; + temp2->isthep = 2; //!< documentation particle + temp2->idhep = id_pair; + for (int j = 0; j < 2; j++) temp2->jmohep[j] = 0; + for (int j = 0; j < 2; j++) temp2->jdahep[j] = 0; + for (int j = 0; j < 5; j++) temp2->phep[j] = 0.0; + temp2->phep[3] += energy; + temp2->phep[4] += mass; + for (int j = 0; j < 4; j++) temp2->vhep[j] = 0.0; + new_event.push_back(*temp2); + + nevhep = read_stdhep(&new_event); + int offset = 2; + + //! read event from lhe file char line[1000]; - - // Find tag bool found_event = false; - while (fgets(line, 1000, lhe_file) != NULL) { - if (strstr(line, " \n", argv[0]); - return 0; - case 'i': - id_beam = atoi(optarg); - break; - case 'j': - id_pair = atoi(optarg); - break; - case '?': - printf("Invalid option or missing option argument; -h to list options\n"); - return 1; - default: - abort(); - } - } - - if (argc - optind < 3) { - printf("Usage: %s [-i beam_id] [-j pair_id] \n", argv[0]); - printf("Use -h for help\n"); - return 1; - } - - int istream = 0; - int ostream = 1; - - // Open input STDHEP file - int n_events = open_read(argv[optind], istream); - if (n_events <= 0) { - fprintf(stderr, "Error: Could not open STDHEP file %s\n", argv[optind]); - return 1; - } - - // Open input LHE file - FILE* lhe_file = fopen(argv[optind + 1], "r"); - if (!lhe_file) { - fprintf(stderr, "Error: Could not open LHE file %s\n", argv[optind + 1]); - close_read(istream); - return 1; + struct stdhep_entry *temp3 = new struct stdhep_entry; + struct stdhep_entry *temp4 = new struct stdhep_entry; + int nup; + fgets(line, 1000, in_file); + sscanf(line, "%d %*d %*f %*f %*f %*f", &nup); + for (int i = 0; i < nup; i++) { + struct stdhep_entry *temp = new struct stdhep_entry; + fgets(line, 1000, in_file); + int icolup0, icolup1; + double phep0 = 100.0; + char blah[1000]; + sscanf(line, "%d %d %d %d %*d %*d %lf %lf %lf %lf %lf %*f %*f", &(temp->idhep), &(temp->isthep), &(temp->jmohep[0]), &(temp->jmohep[1]), &(temp->phep[0]), &(temp->phep[1]), &(temp->phep[2]), &(temp->phep[3]), &(temp->phep[4])); + + // nup == 6 && temp->idhep == 11 for electron/pair of rad + // nup == 7 && temp->idhep == 611 for electron/pair of ap + if((nup == 6 && temp->idhep == 11) || (nup == 7 && temp->idhep == 611)) temp3 = temp; + // nup == 6 && temp->idhep == -11 for positron/pair of rad + // nup == 7 && temp->idhep == -611 for positron/pair of ap + else if((nup == 6 && temp->idhep == -11) || (nup == 7 && temp->idhep == -611)) temp4 = temp; + else delete temp; } - - // Open output STDHEP file - open_write(argv[optind + 2], ostream, n_events); - - printf("Processing events with:\n"); - printf(" Beam particle ID: %d\n", id_beam); - printf(" Pair particle ID: %d\n", id_pair); - - int event_count = 0; - int aprime_count = 0; - int radiative_count = 0; - - // Process each event - while (true) { - // Read STDHEP event - if (!read_next(istream)) { - break; // No more events - } - - vector input_event; - int nevhep = read_stdhep(&input_event); - - // Read LHE event - EventType event_type; - int nup; - LHEParticle scattered_particle, pair_or_electron, positron; - - if (!read_lhe_event(lhe_file, event_type, nup, scattered_particle, - pair_or_electron, positron, id_pair)) { - fprintf(stderr, "Error reading LHE event %d\n", event_count); - break; - } - - // Build output event with unified logic - vector output_event; - - if (event_type == EVENT_APRIME) { - aprime_count++; - - // ===== Find particles in STDHEP by mother relationship ===== - - // Find scattered electron: jmohep[0] = 0 AND idhep = 11 - int scattered_idx = -1; - for (size_t i = 0; i < input_event.size(); i++) { - if (input_event[i].jmohep[0] == 0 && input_event[i].idhep == 11) { - scattered_idx = i; - break; - } - } - - // Find e+e- pair: jmohep[0] = 1 AND idhep = ±11 - int electron_idx = -1; - int positron_idx = -1; - int pair_count = 0; - - for (size_t i = 0; i < input_event.size(); i++) { - if (input_event[i].jmohep[0] == 1) { - if (input_event[i].idhep == 11) { - if (electron_idx == -1) { - electron_idx = i; - pair_count++; - } - } else if (input_event[i].idhep == -11) { - if (positron_idx == -1) { - positron_idx = i; - pair_count++; - } - } - } - } - - // Validation - if (electron_idx == -1 || positron_idx == -1) { - fprintf(stderr, "Error: Event %d - Could not find e+e- pair with jmohep[0]=1\n", event_count); - fprintf(stderr, " Found electron: %d, Found positron: %d\n", - (electron_idx != -1), (positron_idx != -1)); - break; - } - - if (pair_count > 2) { - fprintf(stderr, "Warning: Event %d - Found %d particles with jmohep[0]=1, using first two\n", - event_count, pair_count); - } - - if (input_event.size() > 3 && scattered_idx != -1) { - printf("Info: A-prime event %d has %lu particles (>3), ignoring extras\n", - event_count, input_event.size()); - } - - // Determine vertex for scattered electron (623) - double scattered_vertex[4]; - if (scattered_idx != -1) { - // Use scattered electron's vertex - for (int j = 0; j < 4; j++) { - scattered_vertex[j] = input_event[scattered_idx].vhep[j]; - } - } else { - // No scattered electron found (N=2 case), use (0, 0, 0.02, 0) - scattered_vertex[0] = 0.0; - scattered_vertex[1] = 0.0; - scattered_vertex[2] = 0.02; - scattered_vertex[3] = 0.0; - } - - // ===== Entry 0: Scattered electron (id_beam) ===== - stdhep_entry entry0; - entry0.isthep = 1; - entry0.idhep = id_beam; - entry0.jmohep[0] = 0; - entry0.jmohep[1] = 0; - entry0.jdahep[0] = 0; - entry0.jdahep[1] = 0; - entry0.phep[0] = scattered_particle.phep[0]; - entry0.phep[1] = scattered_particle.phep[1]; - entry0.phep[2] = scattered_particle.phep[2]; - entry0.phep[3] = scattered_particle.phep[3]; - entry0.phep[4] = ELECTRON_MASS; - for (int j = 0; j < 4; j++) { - entry0.vhep[j] = scattered_vertex[j]; - } - output_event.push_back(entry0); - - // ===== Entry 1: A-prime (id_pair) ===== - stdhep_entry entry1; - entry1.isthep = 2; - entry1.idhep = id_pair; - entry1.jmohep[0] = 0; - entry1.jmohep[1] = 0; - entry1.jdahep[0] = 2; - entry1.jdahep[1] = 3; - entry1.phep[0] = pair_or_electron.phep[0]; - entry1.phep[1] = pair_or_electron.phep[1]; - entry1.phep[2] = pair_or_electron.phep[2]; - entry1.phep[3] = pair_or_electron.phep[3]; - entry1.phep[4] = pair_or_electron.phep[4]; - // Use e+e- vertex (displaced vertex) - for (int j = 0; j < 4; j++) { - entry1.vhep[j] = input_event[electron_idx].vhep[j]; - } - output_event.push_back(entry1); - - // ===== Entry 2: Positron (-11) ===== - stdhep_entry entry2 = input_event[positron_idx]; - entry2.jmohep[0] = 2; // Mother is entry 1 (the 622) - entry2.jmohep[1] = 0; - entry2.phep[4] = ELECTRON_MASS; - output_event.push_back(entry2); - - // ===== Entry 3: Electron (11) ===== - stdhep_entry entry3 = input_event[electron_idx]; - entry3.jmohep[0] = 2; // Mother is entry 1 (the 622) - entry3.jmohep[1] = 0; - entry3.phep[4] = ELECTRON_MASS; - output_event.push_back(entry3); - - } else if (event_type == EVENT_RADIATIVE) { - radiative_count++; - - // ===== Find particles in STDHEP by mother relationship ===== - - // Find scattered lepton: jmohep[0] = 0 AND idhep = 13 - int scattered_idx = -1; - for (size_t i = 0; i < input_event.size(); i++) { - if (input_event[i].jmohep[0] == 0 && input_event[i].idhep == 13) { - scattered_idx = i; - break; - } - } - - // Find e+e- pair: jmohep[0] = 1 AND idhep = ±11 - // For radiative, if no particles with jmohep[0]=1, check jmohep[0]=0 - int electron_idx = -1; - int positron_idx = -1; - int pair_count = 0; - - // First try jmohep[0] = 1 - for (size_t i = 0; i < input_event.size(); i++) { - if (input_event[i].jmohep[0] == 1) { - if (input_event[i].idhep == 11) { - if (electron_idx == -1) { - electron_idx = i; - pair_count++; - } - } else if (input_event[i].idhep == -11) { - if (positron_idx == -1) { - positron_idx = i; - pair_count++; - } - } - } - } - - // If not found with jmohep[0]=1, try jmohep[0]=0 (radiative default case) - if (electron_idx == -1 || positron_idx == -1) { - for (size_t i = 0; i < input_event.size(); i++) { - if (input_event[i].jmohep[0] == 0) { - if (input_event[i].idhep == 11 && electron_idx == -1) { - electron_idx = i; - } else if (input_event[i].idhep == -11 && positron_idx == -1) { - positron_idx = i; - } - } - } - } - - // Validation - if (electron_idx == -1 || positron_idx == -1) { - fprintf(stderr, "Error: Event %d - Could not find e+e- pair in radiative event\n", event_count); - fprintf(stderr, " Found electron: %d, Found positron: %d\n", - (electron_idx != -1), (positron_idx != -1)); - break; - } - - if (pair_count > 2) { - fprintf(stderr, "Warning: Event %d - Found %d particles with jmohep[0]=1, using first two\n", - event_count, pair_count); - } - - if (input_event.size() > 2 && scattered_idx != -1) { - printf("Info: Radiative event %d has %lu particles (>2), ignoring extras\n", - event_count, input_event.size()); - } - - // Determine vertex for scattered lepton (623) - double scattered_vertex[4]; - if (scattered_idx != -1) { - // Use scattered lepton's vertex - for (int j = 0; j < 4; j++) { - scattered_vertex[j] = input_event[scattered_idx].vhep[j]; - } - } else { - // No scattered lepton found, use e+e- vertex - for (int j = 0; j < 4; j++) { - scattered_vertex[j] = input_event[electron_idx].vhep[j]; - } - } - - // ===== Entry 0: Scattered lepton (id_beam) ===== - stdhep_entry entry0; - entry0.isthep = 1; - entry0.idhep = id_beam; - entry0.jmohep[0] = 0; - entry0.jmohep[1] = 0; - entry0.jdahep[0] = 0; - entry0.jdahep[1] = 0; - entry0.phep[0] = scattered_particle.phep[0]; - entry0.phep[1] = scattered_particle.phep[1]; - entry0.phep[2] = scattered_particle.phep[2]; - entry0.phep[3] = scattered_particle.phep[3]; - entry0.phep[4] = ELECTRON_MASS; - for (int j = 0; j < 4; j++) { - entry0.vhep[j] = scattered_vertex[j]; - } - output_event.push_back(entry0); - - // ===== Entry 1: Virtual photon (id_pair) ===== - // Calculate 4-momentum sum from STDHEP e+e- (more accurate than LHE) - stdhep_entry entry1; - entry1.isthep = 2; - entry1.idhep = id_pair; - entry1.jmohep[0] = 0; - entry1.jmohep[1] = 0; - entry1.jdahep[0] = 2; - entry1.jdahep[1] = 3; - - // Sum 4-momenta from STDHEP particles - entry1.phep[0] = input_event[electron_idx].phep[0] + input_event[positron_idx].phep[0]; // px - entry1.phep[1] = input_event[electron_idx].phep[1] + input_event[positron_idx].phep[1]; // py - entry1.phep[2] = input_event[electron_idx].phep[2] + input_event[positron_idx].phep[2]; // pz - entry1.phep[3] = input_event[electron_idx].phep[3] + input_event[positron_idx].phep[3]; // E - - // Calculate invariant mass - double E2 = entry1.phep[3] * entry1.phep[3]; - double px2 = entry1.phep[0] * entry1.phep[0]; - double py2 = entry1.phep[1] * entry1.phep[1]; - double pz2 = entry1.phep[2] * entry1.phep[2]; - entry1.phep[4] = sqrt(E2 - px2 - py2 - pz2); - - // Vertex same as e+e- pair - for (int j = 0; j < 4; j++) { - entry1.vhep[j] = input_event[electron_idx].vhep[j]; - } - output_event.push_back(entry1); - - // ===== Entry 2: Electron (11) ===== - stdhep_entry entry2 = input_event[electron_idx]; - entry2.jmohep[0] = 2; // Mother is entry 1 (the 622) - entry2.jmohep[1] = 0; - // Mass already correct (0.000511) - output_event.push_back(entry2); - - // ===== Entry 3: Positron (-11) ===== - stdhep_entry entry3 = input_event[positron_idx]; - entry3.jmohep[0] = 2; // Mother is entry 1 (the 622) - entry3.jmohep[1] = 0; - // Mass already correct (0.000511) - output_event.push_back(entry3); - - } else { - fprintf(stderr, "Error: Unknown event type for event %d\n", event_count); - break; - } - - // Write output event - write_stdhep(&output_event, event_count); - write_file(ostream); - - event_count++; - - if (event_count % 1000 == 0) { - printf("Processed %d events (%d A-prime, %d radiative)\n", - event_count, aprime_count, radiative_count); + + //! Building truth for mother particles + for (int i = 2; i < new_event.size(); i++) { + if (new_event[i].jmohep[0] == 1 + offset) { + new_event[i].jmohep[0] = 2; + if (new_event[1].jdahep[0] == 0) { + new_event[1].jdahep[0] = i+1; + for (int j = 0; j < 4; j++) new_event[1].vhep[j] = new_event[i].vhep[j]; } + new_event[1].phep[0] = temp3->phep[0] + temp4->phep[0]; + new_event[1].phep[1] = temp3->phep[1] + temp4->phep[1]; + new_event[1].phep[2] = temp3->phep[2] + temp4->phep[2]; + new_event[1].phep[3] = temp3->phep[3] + temp4->phep[3]; + new_event[1].phep[4] = sqrt(pow(new_event[1].phep[3], 2) - pow(new_event[1].phep[2], 2) - pow(new_event[1].phep[1], 2) - pow(new_event[1].phep[0], 2)); + + new_event[1].jdahep[1] = i+1; + } + else if (new_event[i].jmohep[0] == 0) { + new_event[i].jmohep[0] = 1; + if (new_event[0].jdahep[0] == 0) new_event[0].jdahep[0] = i+1; + new_event[0].jdahep[1] = i+1; + } } - - // Clean up - close_read(istream); - fclose(lhe_file); - close_write(ostream); - - printf("\nSuccessfully processed %d events\n", event_count); - printf(" A-prime events: %d\n", aprime_count); - printf(" Radiative events: %d\n", radiative_count); - - return 0; + + write_stdhep(&new_event, nevhep); + write_file(ostream); + nevhep++; + + delete temp1; + delete temp2; + delete temp3; + delete temp4; + } } +