From 4d7941d628fcb8d40b10bc1575ea0484f58cc703 Mon Sep 17 00:00:00 2001 From: swenzel Date: Mon, 10 Jun 2024 16:12:20 +0200 Subject: [PATCH 1/2] Standalone event gen executable * easier debugging / valgrinding * mainly to debug the ROOT performance issue https://github.com/root-project/root/issues/15579 --- Detectors/Base/src/Stack.cxx | 3 + run/CMakeLists.txt | 6 ++ run/standalone_eventgen.cxx | 143 +++++++++++++++++++++++++++++++++++ 3 files changed, 152 insertions(+) create mode 100644 run/standalone_eventgen.cxx diff --git a/Detectors/Base/src/Stack.cxx b/Detectors/Base/src/Stack.cxx index de69c866e7b82..aaedd78b13708 100644 --- a/Detectors/Base/src/Stack.cxx +++ b/Detectors/Base/src/Stack.cxx @@ -606,6 +606,9 @@ void Stack::Reset() mTrackRefs->clear(); mTrackIDtoParticlesEntry.clear(); mHitCounter = 0; + mIndexMap.clear(); + mIndexOfPrimaries.clear(); + mTransportedIDs.clear(); } void Stack::Register() diff --git a/run/CMakeLists.txt b/run/CMakeLists.txt index b1d09db4d02c8..7c54c8d443332 100644 --- a/run/CMakeLists.txt +++ b/run/CMakeLists.txt @@ -120,6 +120,12 @@ o2_add_executable(dpl-eventgen SOURCES dpl_eventgen.cxx PUBLIC_LINK_LIBRARIES O2::Generators O2::Framework O2::SimulationDataFormat) +# standalone event generation (no DPL ... mainly for debugging) +o2_add_executable(standalone-eventgen + COMPONENT_NAME sim + SOURCES standalone_eventgen.cxx + PUBLIC_LINK_LIBRARIES O2::Generators O2::SimulationDataFormat) + o2_add_executable(hepmc-publisher COMPONENT_NAME sim SOURCES o2sim_hepmc_publisher.cxx diff --git a/run/standalone_eventgen.cxx b/run/standalone_eventgen.cxx new file mode 100644 index 0000000000000..235541511e8ef --- /dev/null +++ b/run/standalone_eventgen.cxx @@ -0,0 +1,143 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#include "SimulationDataFormat/MCTrack.h" +#include +#include +#include +#include // simple timer from ROOT +#include +#include +#include +#include + +// A simple, non-DPL, executable for event generation + +struct GeneratorTask { + // For readability to indicate where counting certain things (such as events or timeframes) should be of the same order of magnitude + typedef uint64_t GenCount; + std::string generator = "external"; //{"generator", "boxgen", "Name of generator"}; + GenCount eventNum = 3; // {"nEvents", 1, "Number of events"}; + std::string trigger = ""; //{"trigger", "", "Trigger type"}; // + // std::string iniFile = "/home/swenzel/alisw/O2DPG/MC/config/ALICE3/ini/pythia8_pp_136tev.ini"; //{"configFile", "", "INI file containing configurable parameters"}; + std::string iniFile = "/home/swenzel/alisw/O2DPG/MC/config/PWGEM/ini/GeneratorEMCocktail.ini"; + std::string params = ""; // {"configKeyValues", "", "configurable params - configuring event generation internals"}; + long seed = 0; + int aggregate = 10; + std::string vtxModeArg = "kDiamondParam"; + int64_t ttl = -1; // "time-limit", -1, "Maximum run time limit in seconds (default no limit)"}; + std::string outputPrefix = "";// {"output", "", "Optional prefix for kinematics files written on disc. If non-empty, files _Kine.root + _MCHeader.root will be created."}; + GenCount nEvents = 0; + GenCount eventCounter = 0; + GenCount tfCounter = 0; + std::unique_ptr outfile{}; + std::unique_ptr outtree{}; + + // a pointer because object should only be constructed in the device (not during DPL workflow setup) + std::unique_ptr genservice; + TStopwatch timer; + + void init() + { + genservice.reset(new o2::eventgen::GeneratorService); + o2::utils::RngHelper::setGRandomSeed(seed); + nEvents = eventNum; + // helper to parse vertex option; returns true if parsing ok, false if failure + o2::conf::VertexMode vtxmode; + if (!(o2::conf::SimConfig::parseVertexModeString(vtxModeArg, vtxmode))) { + LOG(error) << "Could not parse vtxMode"; + } + + // update config key params + o2::conf::ConfigurableParam::updateFromFile(iniFile); + o2::conf::ConfigurableParam::updateFromString((std::string)params); + // initialize the service + if (vtxmode == o2::conf::VertexMode::kDiamondParam) { + genservice->initService(generator, trigger, o2::eventgen::DiamondParamVertexOption()); + } else if (vtxmode == o2::conf::VertexMode::kNoVertex) { + genservice->initService(generator, trigger, o2::eventgen::NoVertexOption()); + } else if (vtxmode == o2::conf::VertexMode::kCCDB) { + LOG(warn) << "Not yet supported. This needs definition of a timestamp and fetching of the MeanVertex CCDB object"; + } + timer.Start(); + + if (outputPrefix.size() > 0 && !outfile.get()) { + auto kineoutfilename = o2::base::NameConf::getMCKinematicsFileName(outputPrefix.c_str()); + outfile.reset(new TFile(kineoutfilename.c_str(), "RECREATE")); + outtree.reset(new TTree("o2sim", "o2sim")); + } + } + + void run() + { + static int i = 0; + LOG(warn) << "timeframe " << i; + i++; + std::vector mctracks; + o2::dataformats::MCEventHeader mcheader; + std::vector accum; + std::vector accumHeader; + auto mctrack_ptr = &mctracks; + if (outfile.get()) { + auto br = o2::base::getOrMakeBranch(*outtree, "MCTrack", &mctrack_ptr); + br->SetAddress(&mctrack_ptr); + } + + auto toDoThisBatch = std::min((GenCount)aggregate, nEvents - eventCounter); + LOG(info) << "Generating " << toDoThisBatch << " events"; + + for (auto i = 0; i < toDoThisBatch; ++i) { + mctracks.clear(); + genservice->generateEvent_MCTracks(mctracks, mcheader); + // pc.outputs().snapshot(Output{"MC", "MCHEADER", 0}, mcheader); + // pc.outputs().snapshot(Output{"MC", "MCTRACKS", 0}, mctracks); + LOG(info) << "generated " << mctracks.size() << " tracks"; + std::copy(mctracks.begin(),mctracks.end(),std::back_inserter(accum)); + accumHeader.push_back(mcheader); + ++eventCounter; + + // LOG(info) << "mctracks container cap " << mctracks.capacity() << " vs size " << mctracks.size(); + if (outfile.get() && outtree.get()) { + outtree->Fill(); + } + } + + LOG(info) << "Size of tracks accum " << accum.size() * sizeof(o2::MCTrack) << " bytes"; + + // report number of TFs injected for the rate limiter to work + ++tfCounter; + bool time_expired = false; + if (ttl > 0) { + timer.Stop(); + time_expired = timer.RealTime() > ttl; + timer.Start(false); + if (time_expired) { + LOG(info) << "TTL expired after " << eventCounter << " events ... sending end-of-stream"; + } + } + if (eventCounter >= nEvents || time_expired) { + // write out data to disc if asked + if (outfile.get()) { + outtree->SetEntries(eventCounter); + outtree->Write(); + outfile->Close(); + } + } + } // end run +}; // end struct + +int main() { + GeneratorTask task; + task.init(); + for (int i=0; i < task.eventNum/task.aggregate + 1; ++i) { + task.run(); + } +} \ No newline at end of file From bc8f7f19b40459a2d415bfa4a2ab7da67e6669dc Mon Sep 17 00:00:00 2001 From: swenzel Date: Thu, 20 Jun 2024 17:07:43 +0200 Subject: [PATCH 2/2] Improvements for CollisionContextTool (WIP) * important step to create collision contexts for all timeframes in one step * needed to have "pregencollcontext" in O2DPG work with embedding --- .../DigitizationContext.h | 12 +- .../simulation/src/DigitizationContext.cxx | 103 +++++++++++++++++- Steer/src/CollisionContextTool.cxx | 76 ++++++++++++- 3 files changed, 182 insertions(+), 9 deletions(-) diff --git a/DataFormats/simulation/include/SimulationDataFormat/DigitizationContext.h b/DataFormats/simulation/include/SimulationDataFormat/DigitizationContext.h index f531cf6e7f870..ac85c23bbddad 100644 --- a/DataFormats/simulation/include/SimulationDataFormat/DigitizationContext.h +++ b/DataFormats/simulation/include/SimulationDataFormat/DigitizationContext.h @@ -110,6 +110,11 @@ class DigitizationContext /// Check collision parts for vertex consistency. bool checkVertexCompatibility(bool verbose = false) const; + /// retrieves collision context for a single timeframe-id (which may be needed by simulation) + /// (Only copies collision context without QED information. This can be added to the result with the fillQED method + /// in a second step. As a pre-condition, one should have called finalizeTimeframeStructure) + DigitizationContext extractSingleTimeframe(int timeframeid, std::vector const& sources_to_offset); + /// function reading the hits from a chain (previously initialized with initSimChains /// The hits pointer will be initialized (what to we do about ownership??) template @@ -125,8 +130,9 @@ class DigitizationContext // apply collision number cuts and potential relabeling of eventID void applyMaxCollisionFilter(long startOrbit, long orbitsPerTF, int maxColl); - // finalize timeframe structure (fixes the indices in mTimeFrameStartIndex) - void finalizeTimeframeStructure(long startOrbit, long orbitsPerTF); + /// finalize timeframe structure (fixes the indices in mTimeFrameStartIndex) + // returns the number of timeframes + int finalizeTimeframeStructure(long startOrbit, long orbitsPerTF); // Sample and fix interaction vertices (according to some distribution). Makes sure that same event ids // have to have same vertex, as well as event ids associated to same collision. @@ -167,7 +173,7 @@ class DigitizationContext // for each collision we may record/fix the interaction vertex (to be used in event generation) std::vector> mInteractionVertices; - // the collision records _with_ QED interleaved; + // the collision records **with** QED interleaved; std::vector mEventRecordsWithQED; std::vector> mEventPartsWithQED; diff --git a/DataFormats/simulation/src/DigitizationContext.cxx b/DataFormats/simulation/src/DigitizationContext.cxx index f3f40d77042a5..a18f23650dd37 100644 --- a/DataFormats/simulation/src/DigitizationContext.cxx +++ b/DataFormats/simulation/src/DigitizationContext.cxx @@ -19,6 +19,7 @@ #include // for iota #include #include +#include using namespace o2::steer; @@ -196,10 +197,52 @@ o2::parameters::GRPObject const& DigitizationContext::getGRP() const void DigitizationContext::saveToFile(std::string_view filename) const { + // checks if the path content of filename exists ... otherwise it is created before creating the ROOT file + auto ensure_path_exists = [](std::string_view filename) { + try { + // Extract the directory path from the filename + std::filesystem::path file_path(filename); + std::filesystem::path dir_path = file_path.parent_path(); + + // Check if the directory path is empty (which means filename was just a name without path) + if (dir_path.empty()) { + // nothing to do + return true; + } + + // Create directories if they do not exist + if (!std::filesystem::exists(dir_path)) { + if (std::filesystem::create_directories(dir_path)) { + // std::cout << "Directories created successfully: " << dir_path.string() << std::endl; + return true; + } else { + std::cerr << "Failed to create directories: " << dir_path.string() << std::endl; + return false; + } + } + return true; + } catch (const std::filesystem::filesystem_error& ex) { + std::cerr << "Filesystem error: " << ex.what() << std::endl; + return false; + } catch (const std::exception& ex) { + std::cerr << "General error: " << ex.what() << std::endl; + return false; + } + }; + + if (!ensure_path_exists(filename)) { + LOG(error) << "Filename contains path component which could not be created"; + return; + } + TFile file(filename.data(), "RECREATE"); - auto cl = TClass::GetClass(typeid(*this)); - file.WriteObjectAny(this, cl, "DigitizationContext"); - file.Close(); + if (file.IsOpen()) { + auto cl = TClass::GetClass(typeid(*this)); + file.WriteObjectAny(this, cl, "DigitizationContext"); + file.Close(); + } else { + LOG(error) << "Could not write to file " << filename.data(); + } } DigitizationContext* DigitizationContext::loadFromFile(std::string_view filename) @@ -391,13 +434,15 @@ void DigitizationContext::applyMaxCollisionFilter(long startOrbit, long orbitsPe mEventParts = newparts; } -void DigitizationContext::finalizeTimeframeStructure(long startOrbit, long orbitsPerTF) +int DigitizationContext::finalizeTimeframeStructure(long startOrbit, long orbitsPerTF) { mTimeFrameStartIndex = getTimeFrameBoundaries(mEventRecords, startOrbit, orbitsPerTF); LOG(info) << "Fixed " << mTimeFrameStartIndex.size() << " timeframes "; for (auto p : mTimeFrameStartIndex) { LOG(info) << p.first << " " << p.second; } + + return mTimeFrameStartIndex.size(); } std::unordered_map DigitizationContext::getCollisionIndicesForSource(int source) const @@ -483,3 +528,53 @@ void DigitizationContext::sampleInteractionVertices(o2::dataformats::MeanVertexO } } } + +DigitizationContext DigitizationContext::extractSingleTimeframe(int timeframeid, std::vector const& sources_to_offset) +{ + DigitizationContext r; // make a return object + if (mTimeFrameStartIndex.size() == 0) { + LOG(error) << "No timeframe structure determined; Returning empty object. Please call ::finalizeTimeframeStructure before calling this function"; + return r; + } + r.mSimPrefixes = mSimPrefixes; + r.mMuBC = mMuBC; + try { + auto startend = mTimeFrameStartIndex.at(timeframeid); + + auto startindex = startend.first; + auto endindex = startend.second; + + std::copy(mEventRecords.begin() + startindex, mEventRecords.begin() + endindex, std::back_inserter(r.mEventRecords)); + std::copy(mEventParts.begin() + startindex, mEventParts.begin() + endindex, std::back_inserter(r.mEventParts)); + std::copy(mInteractionVertices.begin() + startindex, mInteractionVertices.begin() + endindex, std::back_inserter(r.mInteractionVertices)); + + // let's assume we want to fix the ids for source = source_id + // Then we find the first index that has this source_id and take the corresponding number + // as offset. Thereafter we subtract this offset from all known event parts. + auto perform_offsetting = [&r](int source_id) { + auto indices_for_source = r.getCollisionIndicesForSource(source_id); + int minvalue = r.mEventParts.size(); + for (auto& p : indices_for_source) { + if (p.first < minvalue) { + minvalue = p.first; + } + } + // now fix them + for (auto& p : indices_for_source) { + auto index_into_mEventParts = p.second; + for (auto& part : r.mEventParts[index_into_mEventParts]) { + if (part.sourceID == source_id) { + part.entryID -= minvalue; + } + } + } + }; + for (auto source_id : sources_to_offset) { + perform_offsetting(source_id); + } + + } catch (std::exception) { + LOG(warn) << "No such timeframe id in collision context. Returing empty object"; + } + return r; +} diff --git a/Steer/src/CollisionContextTool.cxx b/Steer/src/CollisionContextTool.cxx index db61343ddf586..5d7fd7d7d345c 100644 --- a/Steer/src/CollisionContextTool.cxx +++ b/Steer/src/CollisionContextTool.cxx @@ -53,6 +53,8 @@ struct Options { bool genVertices = false; // whether to assign vertices to collisions std::string configKeyValues = ""; // string to init config key values long timestamp = -1; // timestamp for CCDB queries + std::string individualTFextraction = ""; // triggers extraction of individuel timeframe components when non-null + // format is path prefix }; enum class InteractionLockMode { @@ -198,7 +200,10 @@ bool parseOptions(int argc, char* argv[], Options& optvalues) "timeframeID", bpo::value(&optvalues.tfid)->default_value(0), "Timeframe id of the first timeframe int this context. Allows to generate contexts for different start orbits")( "first-orbit", bpo::value(&optvalues.firstOrbit)->default_value(0), "First orbit in the run (HBFUtils.firstOrbit)")( "maxCollsPerTF", bpo::value(&optvalues.maxCollsPerTF)->default_value(-1), "Maximal number of MC collisions to put into one timeframe. By default no constraint.")( - "noEmptyTF", bpo::bool_switch(&optvalues.noEmptyTF), "Enforce to have at least one collision")("configKeyValues", bpo::value(&optvalues.configKeyValues)->default_value(""), "Semicolon separated key=value strings (e.g.: 'TPC.gasDensity=1;...')")("with-vertices", "Assign vertices to collisions.")("timestamp", bpo::value(&optvalues.timestamp)->default_value(-1L), "Timestamp for CCDB queries / anchoring"); + "noEmptyTF", bpo::bool_switch(&optvalues.noEmptyTF), "Enforce to have at least one collision")( + "configKeyValues", bpo::value(&optvalues.configKeyValues)->default_value(""), "Semicolon separated key=value strings (e.g.: 'TPC.gasDensity=1;...')")("with-vertices", "Assign vertices to collisions.")("timestamp", bpo::value(&optvalues.timestamp)->default_value(-1L), "Timestamp for CCDB queries / anchoring")( + "extract-per-timeframe", bpo::value(&optvalues.individualTFextraction)->default_value(""), + "Extract individual timeframe contexts. Format required: time_frame_prefix[:comma_separated_list_of_signals_to_offset]"); options.add_options()("help,h", "Produce help message."); @@ -431,7 +436,7 @@ int main(int argc, char* argv[]) // apply max collision per timeframe filters + reindexing of event id (linearisation and compactification) digicontext.applyMaxCollisionFilter(options.tfid * options.orbitsPerTF, options.orbitsPerTF, options.maxCollsPerTF); - digicontext.finalizeTimeframeStructure(options.tfid * options.orbitsPerTF, options.orbitsPerTF); + auto numTimeFrames = digicontext.finalizeTimeframeStructure(options.tfid * options.orbitsPerTF, options.orbitsPerTF); if (options.genVertices) { // TODO: offer option taking meanVertex directly from CCDB ! "GLO/Calib/MeanVertex" @@ -456,5 +461,72 @@ int main(int argc, char* argv[]) } digicontext.saveToFile(options.outfilename); + // extract individual timeframes + if (options.individualTFextraction.size() > 0) { + // we are asked to extract individual timeframe components + + LOG(info) << "Extracting individual timeframe collision contexts"; + // extract prefix path to store these collision contexts + // Function to check the pattern and extract tokens from b + auto check_and_extract_tokens = [](const std::string& input, std::vector& tokens) { + // the regular expression pattern for expected input format + const std::regex pattern(R"(^([a-zA-Z0-9]+)(:([a-zA-Z0-9]+(,[a-zA-Z0-9]+)*))?$)"); + std::smatch matches; + + // Check if the input matches the pattern + if (std::regex_match(input, matches, pattern)) { + // Clear any existing tokens in the vector + tokens.clear(); + + // matches[1] contains the part before the colon which we save first + tokens.push_back(matches[1].str()); + // matches[2] contains the comma-separated list + std::string b = matches[2].str(); + std::regex token_pattern(R"([a-zA-Z0-9]+)"); + auto tokens_begin = std::sregex_iterator(b.begin(), b.end(), token_pattern); + auto tokens_end = std::sregex_iterator(); + + // Iterate over the tokens and add them to the vector + for (std::sregex_iterator i = tokens_begin; i != tokens_end; ++i) { + tokens.push_back((*i).str()); + } + return true; + } + LOG(error) << "Argument for --extract-per-timeframe does not match specification"; + return false; + }; + + std::vector tokens; + if (check_and_extract_tokens(options.individualTFextraction, tokens)) { + auto path_prefix = tokens[0]; + std::vector sources_to_offset{}; + + LOG(info) << "PREFIX is " << path_prefix; + + for (int i = 1; i < tokens.size(); ++i) { + LOG(info) << "Offsetting " << tokens[i]; + sources_to_offset.push_back(digicontext.findSimPrefix(tokens[i])); + } + + // now we are ready to loop over all timeframes + for (int tf_id = 0; tf_id < numTimeFrames; ++tf_id) { + auto copy = digicontext.extractSingleTimeframe(tf_id, sources_to_offset); + + // each individual case gets QED interactions injected + // This should probably be done inside the extraction itself + if (digicontext.isQEDProvided()) { + auto qedSpec = parseInteractionSpec(options.qedInteraction, ispecs, options.useexistingkinematics); + copy.fillQED(qedSpec.name, qedSpec.mcnumberasked, qedSpec.interactionRate); + } + + std::stringstream str; + str << path_prefix << (tf_id + 1) << "/collisioncontext.root"; + copy.saveToFile(str.str()); + LOG(info) << "----"; + copy.printCollisionSummary(options.qedInteraction.size() > 0); + } + } + } + return 0; }