diff --git a/Generators/CMakeLists.txt b/Generators/CMakeLists.txt index f8551a9e5e8ba..7112cada22686 100644 --- a/Generators/CMakeLists.txt +++ b/Generators/CMakeLists.txt @@ -50,6 +50,7 @@ o2_add_library(Generators $<$:src/DecayerPythia8Param.cxx> $<$:src/GeneratorHepMC.cxx> $<$:src/GeneratorHepMCParam.cxx> + $<$:src/AODToHepMC.cxx> PUBLIC_LINK_LIBRARIES FairRoot::Base O2::SimConfig O2::CommonUtils O2::DetectorsBase O2::ZDCBase O2::SimulationDataFormat ${pythia6Target} ${pythiaTarget} ${hepmcTarget} FairRoot::Gen @@ -106,6 +107,7 @@ endif() if(HepMC3_FOUND) list(APPEND headers include/Generators/GeneratorHepMC.h) list(APPEND headers include/Generators/GeneratorHepMCParam.h) + list(APPEND headers include/Generators/AODToHepMC.h) endif() o2_target_root_dictionary(Generators HEADERS ${headers}) diff --git a/Generators/include/Generators/AODToHepMC.h b/Generators/include/Generators/AODToHepMC.h new file mode 100644 index 0000000000000..c3e16f694f746 --- /dev/null +++ b/Generators/include/Generators/AODToHepMC.h @@ -0,0 +1,673 @@ +// Copyright 2023-2099 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. + +/// @author Christian Holm Christensen + +#ifndef ALICEO2_EVENTGEN_AODTOHEPMC_H_ +#define ALICEO2_EVENTGEN_AODTOHEPMC_H_ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace o2 +{ +namespace eventgen +{ +/** + * Convert AOD tables of MC information into a HepMC event structure. + * + * The event structure is kept in memory. + * + * The conventions used here are + * + * - A @e track is a kinematics particle stored in the @c + * o2::aod::McParticles table (aliased as AODToHepMC::Tracks), and + * correspond to a track used during simulation transport. These + * are of type @c o2::aod::McParticle aliased to AODToHepMC::Track. + * + * - A @e generated track is a particle originally made by the event + * generator. These can be recognised by @c + * Track::producedByGenerator() + * + * - A @e particle is a particle stored in the HepMC event structure + * of the class @c HepMC3::GenParticle. + * + * - A @e vertex is where a particle is either produced or disappears. + * A vertex with incoming particles will have a number of outgoing + * particles. Thus a vertex is the @a production vertex for + * out-going particles of that vertex, and the @e end vertex of + * incoming particles of that vertex. Vertexes are of the type @c + * HepMC3::GenVertex. + * + * - The relationship between mother and daughter tracks in + * AODToHepMC::Tracks is encoded by indexes. + * + * - At most two mothers can be stored per track. If a track has + * no mothers, then both indexes are -1. If a track has a + * single mother track, then both indexes are the same. + * + * - A track can have any number of daughters. The first daughter + * track is identified by the first index and the last daughter + * by the second stored index. All tracks indexes with in the + * range from the first to the second index (both inclusive) are + * daughter tracks. If a track has no daughters, then both + * indexes are -1. If a track has one daughter, then both + * indexes are equal. The number of daughters can be obtained + * by + * + * last - first + 1 + * + * if both last and first are zero or greater. + * + * - An event header is the information stored in a row of @c + * o2::aod::McCollisions table (aliased to AODToHepMC::Headers). + * Each event has one such header (aliased as AODToHepMC::Header). + * + * The header stores information on the event such as event number, + * weight, primary interaction point (IP), and so on. + * + * In addition, auxiliary information (impact parameter, + * @f$N_{\mathrm{part}}@f$, and so on) may be stored separate tables. + * + * - The table @c o2::aod::HepMCXSections (aliased to + * AODToHepMC::XSections) stores the event cross section and weight. + * + * - The table @c o2::aod::HepMCPdfInfos (aliased to + * AODToHepMC::PdfInfos) stores the parton distribution function + * parameters used in the event. + * + * - The table @c o2::aod::HepMCHeavyIons (aliased to + * AODToHepMC::PdfHeavyIons) stores the heavy-ion auxiliary + * information, such as the event impact parameter @f$b@f$, + * @f$N_{\mathrm{part}}@f$, @f$N_{\mathrm{coll}}@f$, and so on. + * + * - A HepMC event has a simple header which contains the event + * number, number of particles, and number of vertexes in the event. + * + * Other event information is stored in specific structures. + * + * - The event cross-section(s) and weight(s) is stored in a @c + * HepMC3::GenCrossSection (aliased to + * AODToHepMC::CrossSectionPtr) object. The information to fill + * into that is read from AODToHepMC::XSections table + * + * - The event parton distribution function parameters are stored in + * a @c HepMC3::GenPdfInfo (aliased to AODToHepMC::PdfInfoPtr) + * object. The information to fill into that is read from + * AODToHepMC::PdfInfos table. + * + * - The event heavy-ion parameters are stored in + * a @c HepMC3::GenHeavyIon (aliased to AODToHepMC::HeavyIonPtr) + * object. The information to fill into that is read from + * AODToHepMC::HeavyIons table. + * + * The conversion is done as follows: + * + * - For all MC tracks, create the correspond particle + * (AODToHepMC::makeParticleRecursive) + * + * - Check if we have already created the corresponding particle + * (AODToHepMC::getParticle). If so, then we go on to the next + * track. + * + * - If we are asked to only keep generated tracks + * (AODToHepMC::mOnlyGen set by option @c --hepmc-only-generated), + * i.e., tracks that were posted by the event generator, and this + * track is not such a track (AODToHepMC::isIgnored), nor forced + * (argument @c force) then return nothing. Note that mothers of + * particles that we keep are @e always made (@c force is set to + * true). + * + * - If we do not have a particle yet for the track (look-up in + * AODToHepMC::mParticles) and it is not excluded, then create a + * particle that corresponds to the track, and store the mapping + * from track to particle (in AODToHepMC::mParticles). + * + * - For all mother tracks of the current track, create the + * corresponding particle, following this algorithm (recursion, + * AODToHepMC::makeParticleRecursive). + * + * - If a mother particle has an end vertex, set that vertex as the + * production vertex of the current particle. + * + * - If no mother particle has an end vertex, and this particle is + * not a beam particle and it does have mothers, then create a + * vertex (AODToHepMC::makeVertex) corresponding to the track + * production position, and add this particle as an outgoing + * particle of that vertex. + * + * - In this case, if some mother does not have an outgoing vertex, + * add that mother as an incoming particle to the created + * vertex. + * + * - If the particle is a beam particle (status=4) then store this + * particle as a beam particle. + * + * - If not a beam particle, and the particle has no mothers, mark + * this particle as an orphan. + * + * - Once we have made particle for all tracks, we flesh-out all + * particles. For all tracks (AODToHepMC::fleshOutParticle) + * + * - If this track is ignored (AODToHepMC::isIgnored), or we have no + * particle (AODToHepMC::getParticle) corresponding to this track, + * go on to the next track. + * + * - Get the end vertex of the particle. If any. Set the candidate + * end vertex to this vertex, whether it exists or not. + * + * - Then for each daughter track of the current track, check if it is + * consistent with the end vertex of the particle. + * + * - Check if the dauther is ignored (AODToHepMC::isIgnored). If so, + * move on to the next daughter track. + * + * - Get the particle corresponding to the daughter track + * (AODToHepMC::getParticle.). If no particle is found, move on + * to the next daughter track. + * + * - Check that the daughter particle has an end vertex. If it + * doesn't, mark it as a @a head-less particle. + * + * - If the production vertex of the daughter doesn't match the end + * vertex of the current particle, or the current candidate end + * vertex, then issue a warning, and move on to the next daughter. + * + * - Update the candidate end vertex to the daughter end vertex. + * + * - After processing all daughter particles, and if we have no end + * vertex for the current particle, then + * + * - if we found no candiate end vertex, and the particle is either a + * beam (status=4) or decayed (status=2) particle, issue a + * warning. + * + * - if we do have a candidate end vertex, set that as the end vertex + * of the current particle. + * + * - If, after this, the current particle does have an end vertex, + * loop over all daughters previsouly as head-less and set their + * production vertex to this end vertex. + * + * - At this point, we should have the particles in a proper HepMC + * event structure. + * + * - During simulation transport, the interaction point vertex (IP) + * may not be at (0,0,0,0). Since some consumers of the the HepMC + * event structure may expect the IP to be at (0,0,0,0), we can + * recenter (AODMCToHepMC::recenter) all vertexes of the event. + * This is governed by the member AODMCToHepMC::mRecenter set by the + * option @c --hepmc-recenter + * + * - We can then fill in remaing information into the HepMC event header. + * + * - The event number and weight is set from event header + * (AODToHepMC::makeHeader). + * + * - The event cross section(s and weight(s) is set from + * AODToHepMC::CrossSections table. + * (AODToHepMC::makeXSection). If no AODToHepMC::CrossSections row + * is passed, then dummy values are filled in. + * + * - The event parton distribution function parameters are set from + * AODToHepMC::PdfInfos table. (AODToHepMC::makePdfInfo). If no + * AODToHepMC::PdfInfos row is passed, then dummy values are + * filled in. + * + * - The event heavy-ion parameters are set from + * AODToHepMC::HeavyIons table. (AODToHepMC::makeHeavyIon). If no + * AODToHepMC::HeavyIons row is passed, then dummy values are + * filled in. + * + * - Once all the above is done, we have a complete HepMC event + * (AODToHepMC::mEvent). This event structure is kept in memory. + * + * - Optionally (option @c --hepmc-dumb @e filename) we may write + * the events to a file (AODToHepMC::mOutput). The event is still + * kept in memory. + * + * - Clients of this class (e.g., o2::pwgmm::RivetWrapper) may access + * the event structure for further processing. + * + * The utility @c o2-aod-mc-to-hepmc will read in AODs and write out a + * HepMC event file (plain ASCII). + * + */ +struct AODToHepMC { + /** + * Group of configurables which will be added to an program that + * uses this class. Note that it is really the specialisation of + * framework::OptionManager that propagates the options + * to the program. + */ + struct : framework::ConfigurableGroup { + /** Option for dumping HepMC event structures to disk. Takes one + * argument - the name of the file to write to. */ + framework::Configurable dump{"hepmc-dump", "", + "Dump HepMC event to output"}; + /** Option for only storing particles from the event generator. + * Note, if a particle is stored down, then its mothers will also + * be stored. */ + framework::Configurable onlyGen{"hepmc-only-generated", false, + "Only export generated"}; + /** Use HepMC's tree parsing for building event structure */ + framework::Configurable useTree{"hepmc-use-tree", false, + "Export as tree"}; + /** Floating point precision used when writing to disk */ + framework::Configurable precision{"hepmc-precision", 8, + "Export precision in dump"}; + /** Recenter event at IP=(0,0,0,0). */ + framework::Configurable recenter{"hepmc-recenter", false, + "Recenter the events at (0,0,0,0)"}; + } configs; + /** + * @{ + * @name The containers we subscribe to + */ + /** Alias of MC collisions table */ + using Headers = o2::aod::McCollisions; + /** Alias of MC collisions table row */ + using Header = o2::aod::McCollision; + /** Alias MC particles (tracks) table */ + using Tracks = o2::aod::StoredMcParticles_001; + /** Alias MC particles (tracks) table row */ + using Track = typename Tracks::iterator; + /** Alias auxiliary MC table of cross-sections */ + using XSections = o2::aod::HepMCXSections; + /** Alias auxiliary MC table of parton distribution function + * parameters */ + using PdfInfos = o2::aod::HepMCPdfInfos; + /** Alias auxiliary MC table of heavy-ion parameters */ + using HeavyIons = o2::aod::HepMCHeavyIons; + /** Alias row of auxiliary MC table of cross-sections */ + using XSection = typename XSections::iterator; + /** Alias row of auxiliary MC table of parton distribution function + * parameters */ + using PdfInfo = typename PdfInfos::iterator; + /** Alias row of auxiliary MC table of heavy-ion parameters */ + using HeavyIon = typename HeavyIons::iterator; + /** @} */ + /** + * @{ + * @name Types from HepMC3 + */ + /** Alias HepMC four-vector */ + using FourVector = HepMC3::FourVector; + /** Alias (smart-)pointer to HepMC particle */ + using ParticlePtr = HepMC3::GenParticlePtr; + /** Alias (smart-)pointer to HepMC vertex */ + using VertexPtr = HepMC3::GenVertexPtr; + /** Alias HepMC eventt structure */ + using Event = HepMC3::GenEvent; + /** Alias (smart-)pointer to HepMC heavy-ion object */ + using HeavyIonPtr = HepMC3::GenHeavyIonPtr; + /** Alias (smart-)pointer to HepMC cross section object */ + using CrossSectionPtr = HepMC3::GenCrossSectionPtr; + /** Alias (smart-)pointer to HepMC parton distribution function + * object */ + using PdfInfoPtr = HepMC3::GenPdfInfoPtr; + /** Type used to map tracks to HepMC particles */ + using ParticleMap = std::map; + /** A container of pointers to particles */ + using ParticleVector = std::vector; + /** A container of pointers to vertexes */ + using VertexVector = std::vector; + /** Alias of HepMC writer class */ + using WriterAscii = HepMC3::WriterAscii; + /** The of pointer to HepMC writer class */ + using WriterAsciiPtr = std::shared_ptr; + /** @} */ + /** + * @{ + * @name HepMC3 objects + */ + /** The result of processing */ + Event mEvent; + /** Pointer to cross section-ion information */ + CrossSectionPtr mCrossSec = nullptr; + /** Pointer to heavy-ion information */ + HeavyIonPtr mIon = nullptr; + /** Pointer to parton distribution function information */ + PdfInfoPtr mPdf = nullptr; + /** @} */ + /** + * @{ + * @name Containers etc. + */ + /** Maps tracks to particles */ + ParticleMap mParticles; //! Cache of particles + /** List of vertexes made */ + VertexVector mVertices; //! Cache of vertices + /** List of beam particles */ + ParticleVector mBeams; //! Cache of beam particles + /** Particles without a mother */ + ParticleVector mOrphans; //! Cache of particles w/o mothers + /** @} */ + /** + * @{ + * @name Options and such + */ + /** Output writer, if enabled */ + WriterAsciiPtr mWriter = nullptr; + /** Current sequential event number */ + int mEventNo = 0; + /** The last bunch crossing identifier */ + int mLastBC = -1; + /** If true, only store particles from the generator */ + bool mOnlyGen = false; + /** If true, use HepMC tree parser */ + bool mUseTree = true; + /** Output stream if enabled */ + std::ofstream* mOutput = nullptr; + /** Precision used on the output stream */ + int mPrecision = 16; + /** If true, recenter IP to (0,0,0,0) */ + bool mRecenter = false; + /** @} */ + + /** + * @{ + * @name Interface member functions + */ + /** + * Initialize the converter. Sets internal parameters based on the + * configurables. + */ + virtual void init(); + /** + * Process the collision header and tracks + * + * @param collision Header information + * @param tracks Particle tracks + */ + virtual void process(Header const& collision, Tracks const& tracks); + /** + * Process collision header and HepMC auxiliary information + * + * @param collision Header information + * @param xsections Cross-section table (possible no rows) + * @param pdfs Parton-distribution function table (possible no rows) + * @param heavyions Heavy ion collision table (possible no rows) + */ + virtual void process(Header const& collision, + XSections const& xsections, + PdfInfos const& pdfs, + HeavyIons const& heavyions); + /** + * End of run - closes output file if enabled. This is called via + * specialisation of o2::framework::OutputManager. + */ + virtual bool postRun() + { + enableDump(""); + return true; + } + /** @} */ + protected: + /** + * @{ + * @name Actual worker member functions + */ + /** + * Generate the final event, including fleshing out the vertexes, + * and so on + * + * @param collision Header information + * @param tracks Particle tracks + */ + virtual void makeEvent(Header const& collision, + Tracks const& tracks); + /** + * Set the various fields in the header of the HepMC3::GenEvent + * object + * + * @param header Header object + */ + virtual void makeHeader(Header const& header); + /** + * Make cross-section information. If no entry in the table, + * then make dummy information + */ + virtual void makeXSection(XSections const& xsection); + /** + * Make parton-distribition function information. If no entry + * in the table, then make dummy information + */ + virtual void makePdfInfo(PdfInfos const& pdf); + /** + * Make heavy-ion collision information. If no header given, + * then fill in other reasonable values + */ + virtual void makeHeavyIon(HeavyIons const& heavyion, + Header const& header); + /** + * This is supposed to make the beam particles from the information + * available. However, it seems like we really don't have enough + * information, so for now we will do nothing here. Perhaps the + * user will be forced to provide that information - either via the + * analyser configurables or somehow from somewhere else. + */ + virtual void makeBeams(Header const& header, const VertexPtr ip) {} + /** + * Make all particles. We loop through the MC particles, and for + * each of them create that particle and any possible mother + * particles (recursively). This allows us to traverse the data + * rather straight-forwardly. + * + * @param tracks The MC tracks + */ + virtual void makeParticles(Tracks const& tracks); + /** + * Get particle corresponding to track @a no from particle cache + * + * @param ref Track reference + * @return Pointer to HepMC3::GenParticle or null + */ + virtual ParticlePtr getParticle(Track const& ref) const; + /** + * Check if we are ignoring this track + * + * @param track Track to check */ + virtual bool isIgnored(Track const& track) const; + /** + * Truely make a particle, and its mother particles if any. We add + * vertexes as needed here too. + * + * Note that this traverses the tree from the bottom up. That is, + * we check if a particle has any mothers, and if it does, we create + * those. + * + * However, this can be a bit problematic, since the Kinematic tree + * (and thus the McParticles table) only allows for at most 2 + * mothers. In the case a vertex has 3 or more incoming particles, + * then some of the intermediate particles will be lost. + * + * We remedy that problem by traversing the tree again, but this + * time from the bottom up - that is, we look for daughters of all + * particles, and if they are not registered with the out-going + * vertex, then we reattch the parent to the incoming vertex of the + * daughters. For this to work, we need to map from track index to + * HepMC3::GenParticle. + * + * @param track Current track + * @param tracks MC tracks + * @param force If true, do make the particle even if it would + * otherwise be skipped. + * + * @return Pointer to created particle. + */ + virtual ParticlePtr makeParticleRecursive(Track const& track, + Tracks const& tracks, + bool force = false); + /** + * Generate a HepMC particle from a track. Note that the job here + * is simply to make the object. The more complicated job of adding + * the track to the tree is done above in makeParticleRecursive + * + * @param track MC track + * @param mst Mother status code - updated on return + * @param force Force generation of particle even if onlyGen + * + * @returns Shared pointer to new HepMC particle + */ + virtual ParticlePtr makeParticle(const Track& track, + Int_t& motherStatus, + bool force) const; + /** + * Generate vertex from production vertex of track + * + * @param track MC track + * + * @returns Shared pointer to new HepMC vertex + */ + virtual VertexPtr makeVertex(const Track& track); + /** + * Visit all tracks, but this time look for daughters and + * attach mothers as incoming particles if not already done. + */ + virtual void fleshOut(Tracks const& tracks); + /** + * Flesh out a single particle + * + * @param track The track for which we should flesh out the + * corresponding particle. + */ + virtual void fleshOutParticle(Track const& track, Tracks const& tracks); + /** + * Recenter event to (0,0,0,0). This will use the vertex + * information from the event header to translate all vertexes in + * the event. + * + * @param header Event header + */ + virtual void recenter(Header const& collision); + /** @} */ + /** + * Open dump output, or close if an empty string was given. + * + * @param dump + */ + void enableDump(const std::string& dump); + +}; /** class Generator **/ + +} // namespace eventgen + +namespace framework +{ +/** + * This specialisation of o2::framework::OutputManager ensures that + * we can call the post-processing routine of o2::eventgen::AODToHepMC + * and thus ensure that the possible HepMC is written to disk. + * + * The O2 framework (via o2::framework::adoptAnalysisTask) inspects + * the members of the passed class (@c T) and creates + * o2::framework::OutputManager callbacks for every member. The + * default template for this does nothing. + * + * Thus, to delegate a call to a member of the analysis task (of class + * @c T), we can specialise the @c o2::framework::OutputManager + * template on the @e member type. We will then effectively have + * call-backs for + * + * - @c appendOutput - when the task is constructed + * - @c prepare - when a new set of data is recieved + * - @c finalize - when a set of data has been processed + * - @c postRun - when the run is over + * + * Concretely, we use the @c postRun to flush the HepMC data file + * to disk. + * + * For this to work, the AODToHepMC object must be a member of the + * "Task" class, e.g., + * + * @code + * struct Task { + * o2::eventgen::AODToHepMC mConverter; + * ... + * }; + * + * WorkflowSpec defineDataProcessing(ConfigContext const& cfg) { + * return WorkflowSpec{adaptAnalysisTask(cfg)}; + * } + * @endcode + */ +template <> +struct OutputManager { + /** Type of the target */ + using Target = eventgen::AODToHepMC; + /** Called when task is constructed */ + static bool appendOutput(std::vector&, Target&, uint32_t) { return true; } + /** Called when new data is received */ + static bool prepare(ProcessingContext&, Target&) { return true; } + /** Called when all data has been received */ + static bool postRun(EndOfStreamContext&, Target& t) { return t.postRun(); } + /** Called when the job finishes */ + static bool finalize(ProcessingContext&, Target& t) { return true; } +}; + +/** + * Spacialisation to pull in configurables from the converter. + * + * Ideally, the converter should simply derive from ConfigurableGroup + * and all should flow automatically, but that doesn't work for some + * reason. + * + * For this to work, the AODToHepMC object must be a member of the + * "Task" class, e.g., + * + * @code + * struct Task { + * o2::eventgen::AODToHepMC mConverter; + * ... + * }; + * + * WorkflowSpec defineDataProcessing(ConfigContext const& cfg) { + * return WorkflowSpec{adaptAnalysisTask(cfg)}; + * } + * @endcode + */ +template <> +struct OptionManager { + /** type of the target */ + using Target = eventgen::AODToHepMC; + /** Called when the task is constructed */ + static bool + appendOption(std::vector& options, + Target& target) + { + OptionManager::appendOption(options, target.configs); + return true; + } + /** Called when options are processed */ + static bool + prepare(o2::framework::InitContext& ic, Target& target) + { + OptionManager::prepare(ic, target.configs); + return true; + } +}; + +} // namespace framework +} // namespace o2 + +#endif /* ALICEO2_EVENTGEN_GENERATOR_H_ */ +// Local Variables: +// mode: C++ +// End: diff --git a/Generators/src/AODToHepMC.cxx b/Generators/src/AODToHepMC.cxx new file mode 100644 index 0000000000000..8d508e24ee1cb --- /dev/null +++ b/Generators/src/AODToHepMC.cxx @@ -0,0 +1,584 @@ +// Copyright 2023-2099 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. + +/// @author Christian Holm Christensen + +#include "Generators/AODToHepMC.h" +#include +namespace o2 +{ +namespace eventgen +{ +// ------------------------------------------------------------------- +void AODToHepMC::init() +{ + + mOnlyGen = configs.onlyGen; + mUseTree = configs.useTree; + mPrecision = configs.precision; + mRecenter = configs.recenter; + enableDump(configs.dump); + + LOG(debug) << "\n" + << "=== o2::rivet::Converter ===\n" + << " Dump to output: " << configs.dump << "\n" + << " Only generated tracks: " << mOnlyGen << "\n" + << " Use tree store: " << mUseTree << "\n" + << " Output precision: " << mPrecision; +} +// ------------------------------------------------------------------- +void AODToHepMC::process(Header const& collision, + Tracks const& tracks) +{ + LOG(debug) << "--- Processing track information"; + mBeams.clear(); + mOrphans.clear(); + mParticles.clear(); + mVertices.clear(); + mEvent.clear(); + + makeHeader(collision); + makeParticles(tracks); + makeEvent(collision, tracks); +} +// ------------------------------------------------------------------- +void AODToHepMC::process(Header const& collision, + XSections const& xsections, + PdfInfos const& pdfs, + HeavyIons const& heavyions) +{ + LOG(debug) << "--- Processing auxiliary information"; + makeXSection(xsections); + makePdfInfo(pdfs); + makeHeavyIon(heavyions, collision); +} +// =================================================================== +void AODToHepMC::makeEvent(Header const& collision, + Tracks const& tracks) +{ + LOG(debug) << "Event # " << mEvent.event_number() << " " + << "(# " << mEventNo << " processed)" + << "\n" + << "# of particles " << mParticles.size() << " " + << "(" << tracks.size() << " input tracks)" + << "\n" + << "# of orphans " << mOrphans.size() << "\n" + << "# of beams " << mBeams.size() << "\n" + << "# of vertexes " << mVertices.size(); + if (mUseTree) { + mEvent.reserve(mParticles.size() + mBeams.size(), mVertices.size()); + mEvent.add_tree(mBeams); + } else { + if (mOrphans.size() > 0) { + if (mBeams.size() > 0) { + LOG(debug) << "Event has " << mOrphans.size() << " orphans " + << "out of " << mParticles.size() << " total, " + << "and explicit beams " << mBeams.size(); + } else { + LOG(warning) << "Event has " << mOrphans.size() << " orphans " + << "out of " << mParticles.size() << " total, " + << "but no beams"; + } + // Get collision IP + FourVector v(collision.posX(), + collision.posY(), + collision.posZ(), + collision.t()); + auto ip = std::make_shared(v); + mEvent.add_vertex(ip); + for (auto p : mOrphans) { + ip->add_particle_out(p); + } + // If we have no beam particles, try to make them + if (mBeams.size() == 0) { + makeBeams(collision, ip); + } + } + } + // Flesh out the tracks based on daughter information. + fleshOut(tracks); + if (mWriter) { + // If we have a writer, then dump event to output file + mWriter->write_event(mEvent); + } +} +// ------------------------------------------------------------------- +void AODToHepMC::makeHeader(Header const& header) +{ + int eventNo = header.bcId(); + if (eventNo == mLastBC) { + eventNo = mEventNo; + } + + // Set the event number + mEvent.set_event_number(eventNo); + mEvent.weights().push_back(header.weight()); + LOG(debug) << "Event # " << mEvent.event_number() + << " (BC: " << header.bcId() + << " serial: " << mEventNo + << " last: " << mLastBC << ") " + << " w/weight " << mEvent.weights().front(); + // Store last seen BC + mLastBC = header.bcId(); + // Increase internal counter of event + mEventNo++; +} +// ------------------------------------------------------------------- +void AODToHepMC::makeXSection(XSections const& xsections) +{ + if (not mCrossSec) { + // If we do not have a cross-sections object, create it + mCrossSec = std::make_shared(); + } + + mEvent.set_cross_section(mCrossSec); + mCrossSec->set_cross_section(0.f, 0.f, 0, 0); + + if (xsections.size() <= 0) { + // If we have no info, skip the rest + return; + } + + XSection xsection = xsections.iteratorAt(0); + mCrossSec->set_cross_section(xsection.xsectGen(), + xsection.xsectErr(), + xsection.accepted(), + xsection.attempted()); +} +// ------------------------------------------------------------------- +void AODToHepMC::makePdfInfo(PdfInfos const& pdfs) +{ + if (not mPdf) { + // If we do not have a Parton Distribution Function object, create it + mPdf = std::make_shared(); + } + + mEvent.set_pdf_info(mPdf); + mPdf->set(0, 0, 0.f, 0.f, 0.f, 0.f, 0.f, 0, 0); + + if (pdfs.size() <= 0) { + // If we have no PDF info, skip the rest + return; + } + + PdfInfo pdf = pdfs.iteratorAt(0); + mPdf->set(pdf.id1(), + pdf.id2(), + pdf.x1(), + pdf.x2(), + pdf.scalePdf(), + pdf.pdf1(), + pdf.pdf2(), + pdf.pdfId1(), + pdf.pdfId2()); +} +// ------------------------------------------------------------------- +void AODToHepMC::makeHeavyIon(HeavyIons const& heavyions, + Header const& header) +{ + if (not mIon) { + // Generate heavy ion element if it doesn't exist + mIon = std::make_shared(); + } + + mEvent.set_heavy_ion(mIon); + mIon->impact_parameter = header.impactParameter(); + mIon->event_plane_angle = 0.f; + mIon->Ncoll_hard = 0; + mIon->Npart_proj = 0; + mIon->Npart_targ = 0; + mIon->Nspec_proj_n = 0; + mIon->Nspec_proj_p = 0; + mIon->Nspec_targ_n = 0; + mIon->Nspec_targ_p = 0; + mIon->Ncoll = 0; + mIon->N_Nwounded_collisions = 0; + mIon->Nwounded_N_collisions = 0; + mIon->Nwounded_Nwounded_collisions = 0; + mIon->sigma_inel_NN = 0.f; + mIon->centrality = 0.f; +#ifndef HEPMC3_NO_DEPRECATED + // Deprecated interface with a single eccentricity + mIon->eccentricity = 0.f; +#else + // Newer interface that stores multiple orders of eccentricities. + mIon->eccentricities[1] = 0.f; +#endif + + if (heavyions.size() <= 0) { + // If we have no heavy-ion information, skip the rest + return; + } + + HeavyIon heavyion = heavyions.iteratorAt(0); + float r = 1; + // We need to calculate the ratio projectile to target participants + // so that we may break up the number of spectators and so on. This + // is because the AOD HepMC3HeavyIons table does not store the + // relevant information directly. + if (heavyion.npartProj() < heavyion.npartTarg() and + heavyion.npartTarg() > 0) { + r = heavyion.npartProj() / heavyion.npartTarg(); + } else if (heavyion.npartTarg() < heavyion.npartProj() and + heavyion.npartProj() > 0) { + r = heavyion.npartTarg() / heavyion.npartProj(); + r = (1 - r); + } + + // Heavy ion parameters. Note that number of projectile/target + // proton/neutrons are set by the ratio calculated above. + mIon->impact_parameter = heavyion.impactParameter(); + mIon->event_plane_angle = heavyion.eventPlaneAngle(); + mIon->Ncoll_hard = heavyion.ncollHard(); + mIon->Npart_proj = heavyion.npartProj(); + mIon->Npart_targ = heavyion.npartTarg(); + mIon->Nspec_proj_n = heavyion.spectatorNeutrons() * (1 - r); + mIon->Nspec_proj_p = heavyion.spectatorProtons() * (1 - r); + mIon->Nspec_targ_n = heavyion.spectatorNeutrons() * r; + mIon->Nspec_targ_p = heavyion.spectatorProtons() * r; + mIon->Ncoll = heavyion.ncoll(); + mIon->N_Nwounded_collisions = heavyion.nNwoundedCollisions(); + mIon->Nwounded_N_collisions = heavyion.nwoundedNCollisions(); + mIon->Nwounded_Nwounded_collisions = heavyion.nwoundedNwoundedCollisions(); + mIon->sigma_inel_NN = heavyion.sigmaInelNN(); + mIon->centrality = heavyion.centrality(); +#ifndef HEPMC3_NO_DEPRECATED + mIon->eccentricity = heavyion.eccentricity(); +#else + mIon->eccentricities[1] = heavyion.eccentricity(); +#endif +} +// ------------------------------------------------------------------- +void AODToHepMC::makeParticles(Tracks const& tracks) +{ + for (auto track : tracks) { + makeParticleRecursive(track, tracks); + } +} +// ------------------------------------------------------------------- +AODToHepMC::ParticlePtr AODToHepMC::getParticle(Track const& ref) const +{ + auto iter = mParticles.find(ref.globalIndex()); + if (iter == mParticles.end()) { + return nullptr; + } + + return iter->second; +} +// ------------------------------------------------------------------- +bool AODToHepMC::isIgnored(Track const& track) const +{ + bool fromEG = track.producedByGenerator(); + if (!fromEG and mOnlyGen) { + LOG(debug) << " Track # " << track.globalIndex() + << " from transport, ignored"; + return true; + } + return false; +} +// ------------------------------------------------------------------- +AODToHepMC::ParticlePtr AODToHepMC::makeParticleRecursive(Track const& track, + Tracks const& tracks, + bool force) +{ + ParticlePtr particle = getParticle(track); + + // Check if we already have the particle, and if so, return it + if (particle) { + return particle; + } + + // Make this particle and store it + int motherStatus = 0; + particle = makeParticle(track, motherStatus, force); + if (not particle) { + return nullptr; + } + + // Store mapping from index to particle + mParticles[track.globalIndex()] = particle; + + // Generate mother particles, recurses down tree + ParticleVector mothers; + VertexPtr vout; + for (auto mtrack : track.mothers_as()) { + auto mother = makeParticleRecursive(mtrack, tracks, true); + // If mother not found, continue + if (not mother) { + continue; + } + + // Overrride mother status based on production mechanism of daughter + if (motherStatus != 0) { + mother->set_status(motherStatus); + } + + mothers.push_back(mother); + // Update the production vertex if not set already + if (not vout) { + vout = mother->end_vertex(); + } + } + + // If we have no out vertex, and the particle isn't a beam + // particle, and we have mother particles, then create the + // out-going vertex. + if (not vout and + particle->status() != 4 and + mothers.size() > 0) { + vout = makeVertex(track); + mVertices.push_back(vout); + + // If mothers do not have any end-vertex, add them to the found + // vertex. + for (auto mother : mothers) { + if (not mother->end_vertex()) { + vout->add_particle_in(mother); + } + } + } + + // If we got a out-going vertex, add this particle to that + if (vout) { + vout->add_particle_out(particle); + } + + // If this is a beam particle, add to them + if (particle->status() == 4) { + mBeams.push_back(particle); + } + // if if there no mothers, and this is not beam, then make + // this an orphan. + else if (mothers.size() <= 0) { + mOrphans.push_back(particle); + } + // return the particle + return particle; +} +// ------------------------------------------------------------------- +AODToHepMC::ParticlePtr AODToHepMC::makeParticle(const Track& track, + Int_t& motherStatus, + bool force) const +{ + Int_t pdg = track.pdgCode(); + int hepMCstatus = track.getHepMCStatusCode(); + int egStatus = track.getGenStatusCode(); + int transport = track.getProcess(); + bool fromEG = track.producedByGenerator(); + + // Do not generate particle if it is not from generator and we are + // asked not to make other particles. Note, if a particle has this + // as one of it's mothers, we will make it despite the flag. + if (not force and mOnlyGen and !fromEG) { + return nullptr; + } + + FourVector p(track.px(), + track.py(), + track.pz(), + track.e()); + + // Possibly update mother status depending on the production + // mechanism of this child. + motherStatus = 0; // 200 + uni; // Mother unknown status code + switch (transport) { + case kPDecay: + motherStatus = 2; + break; // Mother decayed! + } + int state = hepMCstatus < 0 ? 200 + transport : hepMCstatus; + + ParticlePtr g = std::make_shared(p, pdg, state); + // g->set_generated_mass(track.GetMass()); + + return g; +} +// ------------------------------------------------------------------- +void AODToHepMC::fleshOut(Tracks const& tracks) +{ + for (auto track : tracks) { + fleshOutParticle(track, tracks); + } +} +// ------------------------------------------------------------------- +void AODToHepMC::fleshOutParticle(Track const& track, Tracks const& tracks) +{ + // If we are only propagating generated tracks, then we need + // not process transported tracks + if (isIgnored(track)) { + return; + } + + // Check that we can find the track in our map + auto particle = getParticle(track); + if (not particle) { + LOG(warning) << "No particle at " << track.globalIndex() << " in map"; + return; + } + + auto endVtx = particle->end_vertex(); + + // If endVtx is null, then the particle has no end vertex. This can + // be because the particle is truly a leaf particle (final-state, + // including after transport), or because the particle wasn't + // properly attached to the event. Since we cannot be sure, in the + // first case, that the status flag accurately reflects the + // sitation, we instead investigate the actual daughters. + // + // We check even if the particle has an end vertex and that all + // daughters have the same production vertex. + int svId = particle->id(); + VertexPtr candidate = endVtx; + ParticleVector headLess; + for (auto dtrack : track.daughters_as()) { + // Check that the daughther is generated by EG. If not, and we + // only store generated tracks, then go on to the next daughter + // (or return?). + if (isIgnored(dtrack)) { + continue; + } + + auto daughter = getParticle(dtrack); + if (not daughter) { + LOG(warning) << "Daughter " << dtrack.globalIndex() + << " of " << track.globalIndex() + << " not found in map!"; + continue; + } + + // We get the production vertex of the daughter. If there's no + // production vertex, then the daughter is deemed "head-less", and + // we will attach it to the end vertex of the mother, if one + // exists or is found below. + auto prodVtx = daughter->production_vertex(); + if (not prodVtx) { + headLess.push_back(daughter); + continue; + } + + // If the mother has an end vertex, but it doesn't match the + // production vertex of the daughter, then this daughter does not + // belong to that mother. This comes about because O2 encodes + // daughters as a range - a la HEPEVT, which requires a specific + // ordering of particles so that all daughters of a vertex are + // consequitive. Since we decide to trust the mother information, + // rather than the daughter information, we will simply disregard + // such daughters. + // + // This check may not be needed + if (endVtx and endVtx->id() != prodVtx->id()) { + continue; + } + + // If we have a current candidate end vertex, but it doesn't match + // the production vertex of the daughter, then we give a warning. + // + // This check may not be needed + if (candidate and prodVtx->id() != candidate->id()) { + LOG(warning) << "Production vertex of daughter " << daughter->id() + << " of " << particle->id() + << " is not the same as previously found from " + << svId; + continue; + } + candidate = prodVtx; + svId = daughter->id(); + } + if (not endVtx) { + // Give warning for decayed or beam particles without and + // end vertex + if (not candidate and + (particle->status() == 4 or + particle->status() == 2)) { + // Only warn for beam and decayed particles + LOG(warning) << "Particle " << track.globalIndex() + << " (" << particle->id() << ")" + << " w/status " << particle->status() + << " does not have an end-vertex, " + << "nor was any found from its daughters"; + } + + // If we have a candidate, set the end vertex + if (candidate) { + endVtx = candidate; + endVtx->add_particle_in(particle); + } + } + + // If we have head-less daughters, add them here + if (endVtx and headLess.size() > 0) { + for (auto daughter : headLess) { + endVtx->add_particle_out(daughter); + } + } +} +// ------------------------------------------------------------------- +void AODToHepMC::enableDump(const std::string& dump) +{ + if (not dump.empty() and mWriter) { + return; + } + + if (dump.empty() and not mWriter) { + return; + } + + if (not dump.empty()) { + LOG(debug) << "o2::rivet::Converter: Open output HepMC file " << dump; + mOutput = new std::ofstream(dump.c_str()); + mWriter = std::make_shared(*mOutput); + mWriter->set_precision(mPrecision); + } else { + LOG(debug) << "o2::river::Converter\n" + << "*********************************\n" + << "Closing output HepMC file\n" + << "*********************************"; + mWriter.reset(); + if (mOutput) { + mOutput->close(); + } + + delete mOutput; + mOutput = nullptr; + } +} +// ------------------------------------------------------------------- +AODToHepMC::VertexPtr AODToHepMC::makeVertex(const Track& track) +{ + FourVector v(track.vx(), + track.vy(), + track.vz(), + track.vt()); + auto vtx = std::make_shared(v); + if (not mUseTree) { + mEvent.add_vertex(vtx); + } + return vtx; +} +// ------------------------------------------------------------------- +void AODToHepMC::recenter(Header const& collision) +{ + FourVector ip(collision.posX(), + collision.posY(), + collision.posZ(), + collision.t()); + + for (auto vertex : mEvent.vertices()) { + vertex->set_position(vertex->position() - ip); + } +} + +// ------------------------------------------------------------------- +} /* namespace eventgen */ +} /* namespace o2 */ +// +// EOF +// diff --git a/run/CMakeLists.txt b/run/CMakeLists.txt index 412c9fddb9086..8c120ee48594d 100644 --- a/run/CMakeLists.txt +++ b/run/CMakeLists.txt @@ -122,6 +122,14 @@ o2_add_executable(mctracks-to-aod SOURCES o2sim_mctracks_to_aod.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::SimulationDataFormat) +o2_add_executable(mc-to-hepmc + COMPONENT_NAME aod + SOURCES o2aod_mc_to_hepmc.cxx + PUBLIC_LINK_LIBRARIES + O2::Framework + O2::SimulationDataFormat + O2::Generators) + o2_add_dpl_workflow(mctracks-to-aod-simple-task SOURCES SimExamples/McTracksToAOD/mctracks_to_aod_simple_task.cxx PUBLIC_LINK_LIBRARIES O2::Framework diff --git a/run/o2aod_mc_to_hepmc.cxx b/run/o2aod_mc_to_hepmc.cxx new file mode 100644 index 0000000000000..091826d8335a7 --- /dev/null +++ b/run/o2aod_mc_to_hepmc.cxx @@ -0,0 +1,238 @@ +// Copyright 2023-2099 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. + +/** @author Christian Holm Christensen */ + +#include +#include +#include +#include + +//-------------------------------------------------------------------- +/** Task to convert AOD MC tables into HepMC event structure + * + * This assumes that the following tables are available on the input: + * + * - @c o2::aod::McCollisions + * - @c o2::aod::McParticles + * - @c o2::aod::HepMCXSections + * - @c o2::aod::HepMCPdfInfos + * - @c o2::aod::HepMCHeavyIons + * + * The application @c o2-sim-mcevent-to-aod publishes these tables. + * + * Ideally, processing auxiliary information should be optional, as + * in @c Task2 below. However, that causes problems. See @c Task2. + */ +struct Task1 { + /** Alias the converter */ + using Converter = o2::eventgen::AODToHepMC; + + /** Our converter */ + Converter mConverter; + + /** @{ + * @name Container types */ + /** Alias converter header table type */ + using Headers = Converter::Headers; + /** Alias converter header type */ + using Header = Converter::Header; + /** Alias converter track table type */ + using Tracks = Converter::Tracks; + /** Alias converter cross-section table type */ + using XSections = Converter::XSections; + /** Alias converter cross-section type */ + using XSection = Converter::XSection; + /** Alias converter parton distribution function table type */ + using PdfInfos = Converter::PdfInfos; + /** Alias converter parton distribution function type */ + using PdfInfo = Converter::PdfInfo; + /** Alias converter heavy-ions table type */ + using HeavyIons = Converter::HeavyIons; + /** Alias converter heavy-ions type */ + using HeavyIon = Converter::HeavyIon; + /** @} */ + + /** Initialize the job */ + void init(o2::framework::InitContext& ic) + { + mConverter.init(); + } + /** Default processing of an event + * + * @param collision Event header + * @param tracks Tracks of the event + */ + void process(Header const& collision, + XSections const& xsections, + PdfInfos const& pdfs, + HeavyIons const& heavyions, + Tracks const& tracks) + { + LOG(debug) << "=== Processing everything ==="; + mConverter.process(collision, + xsections, + pdfs, + heavyions); + mConverter.process(collision, tracks); + } +}; + +//-------------------------------------------------------------------- +/** + * Ideally, this application should work with the case where only + * + * - @c o2::aod::McCollisions + * - @c o2::aod::McParticles + * + * is available, through the use of @c + * o2::framework::ProcessConfigurable, but that seems to fail + * consistently. The issue seems that the application @c + * o2-sim-mcevent-to-aod @c SIGSEGV since it stops publishing the + * tables when the main process of the client (this application) does + * not require those tables. + * + * I tried various combinations of options for @c + * o2-sim-mcevent-to-aod but nothing seems to work. + * + * The error is + * + * @verbatim + * Exception caught: Unable to find OutputSpec with label HepMCXSections. Available Routes: + * - McCollisions: AOD/MCCOLLISION/0 + * - McParticles: AOD/MCPARTICLE/1 + * - : TFF/TFFilename/0 + * - : TFN/TFNumber/0 + * @endverbatim + * + * Interstingly, the application @c o2-sim-mcevent-to-aod works fine + * on its own, e.g., like + * + * @verbatim + * ./o2-sim-kine-publisher \ + * --aggregate-timeframe 1 \ + * --kineFileName pythia8pp | + * ./o2-sim-mcevent-to-aod \ + * --aod-writer-keep dangling + * @endverbatim + * + * works fine. + */ +struct Task2 { + /** Alias the converter type */ + using Converter = o2::eventgen::AODToHepMC; + + /** Our converter */ + Converter mConverter; + + /** @{ + * @name Container types */ + /** Alias converter header table type */ + using Headers = Converter::Headers; + /** Alias converter header type */ + using Header = Converter::Header; + /** Alias converter track table type */ + using Tracks = Converter::Tracks; + /** Alias converter cross-section table type */ + using XSections = Converter::XSections; + /** Alias converter cross-section type */ + using XSection = Converter::XSection; + /** Alias converter parton distribution function table type */ + using PdfInfos = Converter::PdfInfos; + /** Alias converter parton distribution function type */ + using PdfInfo = Converter::PdfInfo; + /** Alias converter heavy-ions table type */ + using HeavyIons = Converter::HeavyIons; + /** Alias converter heavy-ions type */ + using HeavyIon = Converter::HeavyIon; + /** @} */ + + /** Initialize the job */ + void init(o2::framework::InitContext& ic) + { + mConverter.init(); + } + /** Process tracks of an event + * + * @param collision Event header + * @param tracks Tracks of the event + */ + void processTracks(Header const& collision, + Tracks const& tracks) + { + LOG(debug) << "=== Processing event tracks ==="; + mConverter.process(collision, tracks); + } + /** Optional processing of event to extract extra HepMC information + * + * @param collision Event header + * @param xsections Cross-section information + * @param heavyions Heavy ion (geometry) information + */ + void processAux(Header const& collision, + XSections const& xsections, + PdfInfos const& pdfs, + HeavyIons const& heavyions) + { + LOG(debug) << "=== Processing event auxiliaries ==="; + mConverter.process(collision, + xsections, + pdfs, + heavyions); + } + /** Default processing of an event + * + * @param collision Event header + * @param tracks Tracks of the event + */ + void process(Header const& collision, + Tracks const& tracks) + { + LOG(debug) << "=== Processing only tracks ==="; + processTracks(collision, tracks); + } + /** + * Make a process option. + * + * Instead of using the provided preprocessor macro, we instantise + * the template directly here. This is so that we can specify the + * command line argument (@c --hepmc-aux) rather than to rely on an + * auto-generated name (would be @ --processAux). + */ + decltype(o2::framework::ProcessConfigurable{&Task2::processAux, + "hepmc-aux", false, + "Process auxilirary " + "information"}) + doAux = o2::framework::ProcessConfigurable{&Task2::processAux, + "hepmc-aux", false, + "Process auxilirary " + "information"}; +}; + +//-------------------------------------------------------------------- +using WorkflowSpec = o2::framework::WorkflowSpec; +using TaskName = o2::framework::TaskName; +using DataProcessorSpec = o2::framework::DataProcessorSpec; +using ConfigContext = o2::framework::ConfigContext; + +/** Entry point of @a o2-sim-mcevent-to-hepmc */ +WorkflowSpec defineDataProcessing(ConfigContext const& cfg) +{ + using o2::framework::adaptAnalysisTask; + + // Task1: One entry: header, tracks, auxiliary + // Task2: Two entry: header, tracks, and auxiliary + return WorkflowSpec{ + adaptAnalysisTask(cfg, TaskName{"o2-aod-mc-to-hepmc"})}; +} +// +// EOF +//