diff --git a/Detectors/AOD/CMakeLists.txt b/Detectors/AOD/CMakeLists.txt index 5ec99773a5dc9..acd703dcc6be7 100644 --- a/Detectors/AOD/CMakeLists.txt +++ b/Detectors/AOD/CMakeLists.txt @@ -41,14 +41,14 @@ o2_add_executable( workflow COMPONENT_NAME aod-producer TARGETVARNAME targetName - SOURCES src/aod-producer-workflow.cxx src/AODProducerWorkflowSpec.cxx + SOURCES src/aod-producer-workflow.cxx src/AODProducerWorkflowSpec.cxx src/AODMcProducerHelpers.cxx PUBLIC_LINK_LIBRARIES internal::AODProducerWorkflow O2::Version ) o2_add_executable( workflow COMPONENT_NAME aod-mc-producer - SOURCES src/aod-mc-producer-workflow.cxx src/AODMcProducerWorkflowSpec.cxx + SOURCES src/aod-mc-producer-workflow.cxx src/AODMcProducerWorkflowSpec.cxx src/AODMcProducerHelpers.cxx PUBLIC_LINK_LIBRARIES internal::AODProducerWorkflow O2::Version ) diff --git a/Detectors/AOD/include/AODProducerWorkflow/AODMcProducerHelpers.h b/Detectors/AOD/include/AODProducerWorkflow/AODMcProducerHelpers.h new file mode 100644 index 0000000000000..42431d19cb210 --- /dev/null +++ b/Detectors/AOD/include/AODProducerWorkflow/AODMcProducerHelpers.h @@ -0,0 +1,324 @@ +// 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. + +/// @file AODMcProducerHelpers.h +/// @author Christian Holm Christensen +/// common helpers for AOD MC producers + +#ifndef O2_AODMCPRODUCER_HELPERS +#define O2_AODMCPRODUCER_HELPERS +#include +#include +#include +#include +#include +#include +#include +#include + +/** + * Utilities to transform simulated data into AO2D tables. + * + * The function templates below are templated on the cursor type over + * the relevant AOD tables. Such a table can be obtained from the + * ProcessingContext @c pc + * + * @code + * auto builder = pc.make(OutputForTable::ref()); + * auto cursor = builder->cursor
(); + * @endcode + * + * If the task uses the @c Produces
template, + * + * @code + * Produces
mTable; + * @endcode + * + * then a cursor is obtained via, + * + * @code + * auto cursor = mTable.cursor; + * @endcode + * + * Note that these functions cannot be moved into a compilation unit, + * because that would require deducing the table cursor type, by + * f.ex. + * + * @code + * template + * struct TableCursor { + * using cursor_t = decltype(std::declval() + * .cursor
()); + * }; + * using CollisionCursor = TableCursor:cursor_t; + * @endcode + * + * but since cursors are really Lambdas and Lambda types are specific + * to the compilation unit, then the implementation file (compilation + * unit) of these functions definitions and their use (another + * compilation unit) would have different types of the the cursers, + * and thus not be able to link. More information is given at + * https://stackoverflow.com/questions/50033797. + */ +namespace o2::aodmchelpers +{ +//================================================================== +/** + * Deduce cursor type and wrap in std::function + */ +template +struct TableCursor { + using type = decltype(framework::FFL(std::declval() + .cursor
())); +}; +//================================================================== +/** Cursor over aod::McCollisions */ +using CollisionCursor = TableCursor::type; +/** Cursor over aod::McParticles */ +using ParticleCursor = TableCursor::type; +/** Cursor over aod::HepMCXSections */ +using XSectionCursor = TableCursor::type; +/** Cursor over aod::HepMCPdfInfos */ +using PdfInfoCursor = TableCursor::type; +/** Cursor over aod::HepMCHeavyIons */ +using HeavyIonCursor = TableCursor::type; +//================================================================== +/** Types of updates on HepMC tables. */ +enum HepMCUpdate { + never, + always, + anyKey, + allKeys +}; + +//================================================================== +/** + * Check if header has keys. If the argument @a anyNotAll is true, + * then this member function returns true if @e any of the keys + * were found. If @a anyNotAll is false, then return true only if + * @a all keys were found. + * + * @param header MC event header + * @param keys Keys to look for + * @param anyNotAll If true, return true if @e any key was found. + * If false, return true only if @a all keys were found + * + * @return true if any or all keys were found + */ +bool hasKeys(o2::dataformats::MCEventHeader const& header, + const std::vector& keys, + bool anyNotall = true); +//-------------------------------------------------------------------- +/** + * Get a property from the header, or if not set or not valid, a + * default value. + * + * @param header The MC event header + * @param key Key to look for + * @param def Value to return if key is not found + * + * @return Value of key or def if key is not found + */ +template +const T getEventInfo(o2::dataformats::MCEventHeader const& header, + std::string const& key, + T const& def) +{ + if (not header.hasInfo(key)) + return def; + + bool isValid = false; + const T& val = header.getInfo(key, isValid); + if (not isValid) + return def; + + return val; +} +//==================================================================== +/** + * Fill in collision information. This is read from the passed MC + * header and stored in the MCCollision table. The member function + * returns the encoded generator ID. + * + * @param cursor Cursor over o2::aod::McCollisions table + * @param bcId Bunch-crossing Identifier + * @param time Time of collisions + * @param header Event header from generator + * @param generatorId Default generator + * @param sourceId Identifier of source + * + * @return encoded generator ID + */ +short updateMCCollisions(const CollisionCursor& cursor, + int bcId, + float time, + o2::dataformats::MCEventHeader const& header, + short generatorId = 0, + int sourceId = 0, + unsigned int mask = 0xFFFFFFF0); +//-------------------------------------------------------------------- +/** + * Fill in HepMC cross-section table from event generator header. + * + * @param cursor Cursor over o2::aod::HepMCXSections table + * @param collisionID Identifier of collision (as given updateMCCollision) + * @param generatorID Encoded generator ID + * @param header Event header from generator + * @param anyNotAll If true, then any key present trigger and update. + * If false, then all keys must be present to update + * the table. + * + * @return true if table was updated + */ +bool updateHepMCXSection(const XSectionCursor& cursor, + int collisionID, + short generatorID, + o2::dataformats::MCEventHeader const& header, + HepMCUpdate when = HepMCUpdate::anyKey); +//-------------------------------------------------------------------- +/** + * Fill in HepMC parton distribution function table from event + * generator header + * + * @param cursor Cursor over o2::aod::HepMCXSections table + * @param collisionID Identifier of collision (as given updateMCCollision) + * @param generatorID Encoded generator ID + * @param header Event header from generator + * @param anyNotAll If true, then any key present trigger and update. + * If false, then all keys must be present to update + * the table. + * + * @return true if table was updated + */ +bool updateHepMCPdfInfo(const PdfInfoCursor& cursor, + int collisionID, + short generatorID, + o2::dataformats::MCEventHeader const& header, + HepMCUpdate when = HepMCUpdate::anyKey); +//-------------------------------------------------------------------- +/** + * Fill in HepMC heavy-ion table from generator event header. + * + * @param cursor Cursor over o2::aod::HepMCXSections table + * @param collisionID Identifier of collision (as given updateMCCollision) + * @param generatorID Encoded generator ID + * @param header Event header from generator + * @param anyNotAll If true, then any key present trigger and update. + * If false, then all keys must be present to update + * the table. + * + * @return true if table was updated + */ +bool updateHepMCHeavyIon(const HeavyIonCursor& cursor, + int collisionID, + short generatorID, + o2::dataformats::MCEventHeader const& header, + HepMCUpdate when = HepMCUpdate::anyKey); +//-------------------------------------------------------------------- +/** + * Type of mapping from track number to row index + */ +using TrackToIndex = std::unordered_map; +//-------------------------------------------------------------------- +/** + * Update aod::McParticles table with information from an MC track. + * + * @param cursor Cursor over aod::McParticles table + * @param mapping Maps track number to index in table + * @param collisionID Collision identifier + * @param track Track to update table with + * @param tracks List of all tracks of current collision + * @param flags Base flags of this track + * @param weightMask Mask on weight floating point value + * @param momentumMask Mask on momentum floating point values + * @param positionMask Mask on position floating point values + */ +void updateParticle(const ParticleCursor& cursor, + const TrackToIndex& toStore, + int collisionID, + o2::MCTrack const& track, + std::vector const& tracks, + uint8_t flags = 0, + uint32_t weightMask = 0xFFFFFFF0, + uint32_t momentumMask = 0xFFFFFFF0, + uint32_t positionMask = 0xFFFFFFF0); +//-------------------------------------------------------------------- +/** + * Update aod::McParticles table with tracks from MC. + * + * To add particles from many events, one will do + * + * @code + * TrackToIndex preselect = findMcTracksToStore(...); + * + * size_t offset = 0; + * for (auto event : events) + * offset = updateParticles(cursor, + * event.getCollisionID(), + * event.getTracks(), + * offset, + * filter, + * event.isBackground(), + * preselect); + * @endcode + * + * Here @a preselect must be a map from track number to a positive + * value. Tracks that are mapped as such in @a preselect are stored + * in addition to other tracks selected by the function. Note that @a + * preselect may be empty. + * + * If @a filter is false, then @a all tracks will be stored. + * + * If @a filter is true, then tracks that are + * + * - generated by the generator, + * - physical primaries + * (MCTrackNavigator::isPhysicalPrimary), + * - to be kept for physics + * (MCTrackNavigator::isKeepPhysics), or + * - is listed with a positive value in @a preselect, or + * - either a mother or daughter of one such track, then + * + * that track is kept + * + * On return, the @a preselect will map from track number (index in + * the @a tracks container) to the table row index (including offset + * from previous events in the same time-frame). + * + * @param cursor Cursor over aod::McParticles + * @param int Collision identifier + * @param tracks List of all tracks of current collision + * @param offset Index just beyond last table entry + * @param filter Filter tracks + * @param background True of from background event + * @param preselect Mapping of preselected tracks + * @param weightMask Mask on weight floating point value + * @param momentumMask Mask on momentum floating point values + * @param positionMask Mask on position floating point values + * + * @return Index beyond the last particle added to table + */ +uint32_t updateParticles(const ParticleCursor& cursor, + int collisionID, + std::vector const& tracks, + TrackToIndex& preselect, + uint32_t offset = 0, + bool filter = false, + bool background = false, + uint32_t weightMask = 0xFFFFFFF0, + uint32_t momentumMask = 0xFFFFFFF0, + uint32_t positionMask = 0xFFFFFFF0); +} // namespace o2::aodmchelpers + +#endif /* O2_AODMCPRODUCER_HELPERS */ +// Local Variables: +// mode: C++ +// End: diff --git a/Detectors/AOD/include/AODProducerWorkflow/AODMcProducerWorkflowSpec.h b/Detectors/AOD/include/AODProducerWorkflow/AODMcProducerWorkflowSpec.h index 392dad61df4b3..674fbf7bfcd05 100644 --- a/Detectors/AOD/include/AODProducerWorkflow/AODMcProducerWorkflowSpec.h +++ b/Detectors/AOD/include/AODProducerWorkflow/AODMcProducerWorkflowSpec.h @@ -15,14 +15,14 @@ #define O2_AODMCPRODUCER_WORKFLOW_SPEC #include "AODProducerHelpers.h" +#include "AODMcProducerHelpers.h" #include "CommonDataFormat/InteractionRecord.h" #include "Framework/AnalysisDataModel.h" #include "Framework/AnalysisHelpers.h" -#include "Framework/DataProcessorSpec.h" #include "Framework/Task.h" +#include "Framework/DataProcessorSpec.h" #include "Steer/MCKinematicsReader.h" -#include "TMap.h" -#include "TStopwatch.h" +#include #include #include @@ -41,30 +41,95 @@ class AODMcProducerWorkflowDPL : public Task void run(ProcessingContext& pc) final; void endOfStream(EndOfStreamContext& ec) final; + /** Some types we will use */ + using MCEventHeader = dataformats::MCEventHeader; + using McCollisions = aod::McCollisions; + using McParticles = aod::StoredMcParticles; + using Origins = aod::Origins; + using XSections = aod::HepMCXSections; + using PdfInfos = aod::HepMCPdfInfos; + using HeavyIons = aod::HepMCHeavyIons; + using MCKinematicsReader = steer::MCKinematicsReader; + /** */ private: + /** some other types we will use */ + using CollisionCursor = aodmchelpers::CollisionCursor; + using XSectionCursor = aodmchelpers::XSectionCursor; + using PdfInfoCursor = aodmchelpers::PdfInfoCursor; + using HeavyIonCursor = aodmchelpers::HeavyIonCursor; + using HepMCUpdate = aodmchelpers::HepMCUpdate; + /** + * Update the header (collision and HepMC aux) information. + * + * When updating the HepMC aux tables, we take the relevant policies + * into account (mXSectionUpdate, mPdfInfoUpdate, mHeavyIonUpdate). + * + * - If a policy is "never", then the corresponding table is never + * updated. + * + * - If the policy is "always", then the table is always + * update. + * + * - If the policy is either "anyKey" or "allKeys", _and_ + * this is the first event, then we check if any or all keys, + * respectively are present in the header. + * + * - If that check fails, then we do not update and set the + * corresponding policy to be "never". + * + * - If the check succeeds, then we do update the table, and set + * the corresponding policty to "always". + * + * In this way, we will let the first event decide what to do for + * subsequent events and thus avoid too many string comparisions. + * + * @param collisionCursor Cursor over aod::McCollisions + * @param xSectionCursor Cursor over aod::HepMCXSections + * @param pdfInfoCursor Cursor over aod::HepMCPdfInfos + * @param heavyIonCursor Cursor over aod::HepMCHeavyIons + * @param header Header to read information from + * @param collisionID Index of collision in table + * @param bcID Current event identifier (bcID) + * @param time Time of event + * @param generatorID Generator identifier, if any + * @param sourceID Source identifier + * + */ + void updateHeader(CollisionCursor& collisionCursor, + XSectionCursor& xSectionCursor, + PdfInfoCursor& pdfInfoCursor, + HeavyIonCursor& heavyIonCursor, + const MCEventHeader& header, + int collisionID, // Index + int bcID, + float time, + short generatorID, + int sourceID); + /** Current timeframe number */ int64_t mTFNumber{1}; - int mTruncate{1}; + /** Whether to filter tracks so that only primary particles, + particles from EG, or their parents are written out */ int mFilterMC{0}; + /** Whether to enable embedding */ bool mEnableEmbed{false}; - std::string mMCHeaderFNames; + /** LPM production tag, for anchoring */ TString mLPMProdTag{""}; + /** Pass identifier for anchoring */ TString mAnchorPass{""}; + /** Production identifier for anchoring */ TString mAnchorProd{""}; + /** Reconstruction pass for anchoring */ TString mRecoPass{""}; + /** Timer */ TStopwatch mTimer; + /** Rules for when to update HepMC tables */ + HepMCUpdate mXSectionUpdate = HepMCUpdate::anyKey; + HepMCUpdate mPdfInfoUpdate = HepMCUpdate::anyKey; + HepMCUpdate mHeavyIonUpdate = HepMCUpdate::anyKey; std::string mSimPrefix; - o2::aodhelpers::TripletsMap_t mToStore; - - // keep track event/source id for each mc-collision - // using map and not unordered_map to ensure - // correct ordering when iterating over container elements - std::vector> mMCColToEvSrc; - // MC production metadata holder - std::vector mMetaDataKeys; - std::vector mMetaDataVals; bool mIsMDSent{false}; // truncation is enabled by default @@ -72,10 +137,6 @@ class AODMcProducerWorkflowDPL : public Task uint32_t mMcParticleW = 0xFFFFFFF0; // 19 bits uint32_t mMcParticlePos = 0xFFFFFFF0; // 19 bits uint32_t mMcParticleMom = 0xFFFFFFF0; // 19 bits - - template - void fillMCParticlesTable(o2::steer::MCKinematicsReader& mcReader, - const MCParticlesCursorType& mcParticlesCursor); }; /// create a processor spec diff --git a/Detectors/AOD/include/AODProducerWorkflow/AODProducerHelpers.h b/Detectors/AOD/include/AODProducerWorkflow/AODProducerHelpers.h index f53b2a3082260..dc6f589977cda 100644 --- a/Detectors/AOD/include/AODProducerWorkflow/AODProducerHelpers.h +++ b/Detectors/AOD/include/AODProducerWorkflow/AODProducerHelpers.h @@ -20,6 +20,7 @@ #include #include #include +#include namespace o2::aodhelpers { @@ -48,6 +49,15 @@ struct TripletEqualTo { typedef boost::unordered_map TripletsMap_t; +template +framework::Produces createTableCursor(framework::ProcessingContext& pc) +{ + framework::Produces c; + c.resetCursor(pc.outputs() + .make(framework::OutputForTable::ref())); + c.setLabel(o2::aod::MetadataTrait::metadata::tableLabel()); + return c; +} } // namespace o2::aodhelpers #endif /* O2_AODPRODUCER_HELPERS */ diff --git a/Detectors/AOD/include/AODProducerWorkflow/AODProducerWorkflowSpec.h b/Detectors/AOD/include/AODProducerWorkflow/AODProducerWorkflowSpec.h index caae2fd6b9268..13a55526afc56 100644 --- a/Detectors/AOD/include/AODProducerWorkflow/AODProducerWorkflowSpec.h +++ b/Detectors/AOD/include/AODProducerWorkflow/AODProducerWorkflowSpec.h @@ -14,6 +14,7 @@ #ifndef O2_AODPRODUCER_WORKFLOW_SPEC #define O2_AODPRODUCER_WORKFLOW_SPEC +#include "AODMcProducerHelpers.h" #include "DataFormatsEMCAL/Cell.h" #include "DataFormatsGlobalTracking/RecoContainer.h" #include "DataFormatsPHOS/Cell.h" @@ -63,11 +64,19 @@ class BunchCrossings /// return the sorted vector of increaing BC times std::vector const& getBCTimeVector() const { return mBCTimeVector; } - /// Performs a "lower bound" search for timestamp within the bunch crossing data. - /// Returns the smallest bunch crossing (index and value) equal or greater than timestamp. - /// The functions is expected to perform much better than - /// a binary search in the bunch crossing data directly. Expect O(1) instead of O(log(N)) at the cost - /// of the additional memory used by this class. + /// Performs a "lower bound" search for timestamp within the bunch + /// crossing data. + /// + /// Returns the smallest bunch crossing (index and value) equal or + /// greater than timestamp. + /// + /// The functions is expected to perform much better than a binary + /// search in the bunch crossing data directly. Expect O(1) instead + /// of O(log(N)) at the cost of the additional memory used by this + /// class. + /// + /// This is _not_ O(1). The loop below makes it at least O(N). The + /// call to std::lower_bound is O(log(N)). std::pair lower_bound(uint64_t timestamp) const { // a) determine the timewindow @@ -218,15 +227,6 @@ class AODProducerWorkflowDPL : public Task return std::uint64_t(mStartIR.toLong()) + relativeTime_to_LocalBC(relativeTimeStampInNS); } - template - Produces createTableCursor(ProcessingContext& pc) - { - Produces c; - c.resetCursor(pc.outputs().make(OutputForTable::ref())); - c.setLabel(o2::aod::MetadataTrait::metadata::tableLabel()); - return c; - } - bool mPropTracks{false}; bool mPropMuons{false}; float mTrackQCFraction{0.00}; @@ -541,9 +541,68 @@ class AODProducerWorkflowDPL : public Task template void fillStrangenessTrackingTables(const o2::globaltracking::RecoContainer& data, V0C& v0Cursor, CC& cascadeCursor, D3BC& decay3bodyCursor); - template + /** some other types we will use */ + using MCCollisionCursor = aodmchelpers::CollisionCursor; + using XSectionCursor = aodmchelpers::XSectionCursor; + using PdfInfoCursor = aodmchelpers::PdfInfoCursor; + using HeavyIonCursor = aodmchelpers::HeavyIonCursor; + using MCParticlesCursor = aodmchelpers::ParticleCursor; + using HepMCUpdate = aodmchelpers::HepMCUpdate; + using MCEventHeader = dataformats::MCEventHeader; + /** Rules for when to update HepMC tables */ + HepMCUpdate mXSectionUpdate = HepMCUpdate::anyKey; + HepMCUpdate mPdfInfoUpdate = HepMCUpdate::anyKey; + HepMCUpdate mHeavyIonUpdate = HepMCUpdate::anyKey; + /** + * Update the header (collision and HepMC aux) information. + * + * When updating the HepMC aux tables, we take the relevant policies + * into account (mXSectionUpdate, mPdfInfoUpdate, mHeavyIonUpdate). + * + * - If a policy is "never", then the corresponding table is never + * updated. + * + * - If the policy is "always", then the table is always + * update. + * + * - If the policy is either "anyKey" or "allKeys", _and_ + * this is the first event, then we check if any or all keys, + * respectively are present in the header. + * + * - If that check fails, then we do not update and set the + * corresponding policy to be "never". + * + * - If the check succeeds, then we do update the table, and set + * the corresponding policty to "always". + * + * In this way, we will let the first event decide what to do for + * subsequent events and thus avoid too many string comparisions. + * + * @param collisionCursor Cursor over aod::McCollisions + * @param xSectionCursor Cursor over aod::HepMCXSections + * @param pdfInfoCursor Cursor over aod::HepMCPdfInfos + * @param heavyIonCursor Cursor over aod::HepMCHeavyIons + * @param header Header to read information from + * @param collisionID Index of collision in the table + * @param bcID Current event identifier (bcID) + * @param time Time of event + * @param generatorID Generator identifier, if any + * @param sourceID Source identifier + * + */ + void updateMCHeader(MCCollisionCursor& collisionCursor, + XSectionCursor& xSectionCursor, + PdfInfoCursor& pdfInfoCursor, + HeavyIonCursor& heavyIonCursor, + const MCEventHeader& header, + int collisionID, + int bcID, + float time, + short generatorID, + int sourceID); + void fillMCParticlesTable(o2::steer::MCKinematicsReader& mcReader, - MCParticlesCursorType& mcParticlesCursor, + MCParticlesCursor& mcParticlesCursor, const gsl::span& primVer2TRefs, const gsl::span& GIndices, const o2::globaltracking::RecoContainer& data, diff --git a/Detectors/AOD/src/AODMcProducerHelpers.cxx b/Detectors/AOD/src/AODMcProducerHelpers.cxx new file mode 100644 index 0000000000000..492ac93e10ce2 --- /dev/null +++ b/Detectors/AOD/src/AODMcProducerHelpers.cxx @@ -0,0 +1,431 @@ +// 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. + +/// @file AODMcProducerHelpers.h +/// @author Christian Holm Christensen +/// common helpers for AOD MC producers +#include "AODProducerWorkflow/AODMcProducerHelpers.h" +#include +#include +#include +#include +#define verbosity debug + +namespace o2::aodmchelpers +{ +//================================================================== +bool hasKeys(o2::dataformats::MCEventHeader const& header, + const std::vector& keys, + bool anyNotAll) +{ + auto check = [&header](const std::string& key) { // Do not format + return header.hasInfo(key); + }; + return (anyNotAll ? // Do not format + std::any_of(keys.cbegin(), keys.cend(), check) + : // Do not format + std::all_of(keys.cbegin(), keys.cend(), check)); +} +//==================================================================== +short updateMCCollisions(const CollisionCursor& cursor, + int bcId, + float time, + o2::dataformats::MCEventHeader const& header, + short generatorId, + int sourceId, + unsigned int mask) +{ + using Key = o2::dataformats::MCInfoKeys; + using GenProp = o2::mcgenid::GeneratorProperty; + using namespace o2::math_utils; + + int subGenId = getEventInfo(header, GenProp::SUBGENERATORID, -1); + int genId = getEventInfo(header, GenProp::GENERATORID, + int(generatorId)); + float weight = getEventInfo(header, Key::weight, 1.f); + + auto encodedGeneratorId = o2::mcgenid::getEncodedGenId(genId, + sourceId, + subGenId); + + LOG(verbosity) << "Updating MC Collisions table w/bcId=" << bcId; + cursor(0, + bcId, + encodedGeneratorId, + truncateFloatFraction(header.GetX(), mask), + truncateFloatFraction(header.GetY(), mask), + truncateFloatFraction(header.GetZ(), mask), + truncateFloatFraction(time, mask), + truncateFloatFraction(weight, mask), + header.GetB()); + return encodedGeneratorId; +} +//-------------------------------------------------------------------- +bool updateHepMCXSection(const XSectionCursor& cursor, + int collisionID, + short generatorID, + o2::dataformats::MCEventHeader const& header, + HepMCUpdate when) +{ + using Key = o2::dataformats::MCInfoKeys; + + if (when == HepMCUpdate::never or + (when != HepMCUpdate::always and + not hasKeys(header, { // Do not + Key::acceptedEvents, // mess with + Key::attemptedEvents, // with + Key::xSection, // my + Key::xSectionError}, // formatting + when == HepMCUpdate::anyKey))) { + return false; + } + + LOG(verbosity) << "Updating HepMC cross-section table w/collisionId " + << collisionID; + cursor(0, + collisionID, + generatorID, + getEventInfo(header, Key::acceptedEvents, 0), + getEventInfo(header, Key::attemptedEvents, 0), + getEventInfo(header, Key::xSection, 0.f), + getEventInfo(header, Key::xSectionError, 0.f), + getEventInfo(header, "ptHard", 1.f), + getEventInfo(header, "MPI", -1), + getEventInfo(header, "processId", -1)); + return true; +} +//-------------------------------------------------------------------- +bool updateHepMCPdfInfo(const PdfInfoCursor& cursor, + int collisionID, + short generatorID, + o2::dataformats::MCEventHeader const& header, + HepMCUpdate when) +{ + using Key = o2::dataformats::MCInfoKeys; + + if (when == HepMCUpdate::never or + (when != HepMCUpdate::always and // Do + not hasKeys(header, {Key::pdfParton1Id, // not + Key::pdfParton2Id, // mess + Key::pdfCode1, // with + Key::pdfCode2, // my + Key::pdfX1, // formatting + Key::pdfX2, // . + Key::pdfScale, // It + Key::pdfXF1, // is + Key::pdfXF2}, // better + when == HepMCUpdate::anyKey))) { + return false; + } + + LOG(verbosity) << "Updating HepMC PDF table w/collisionId " << collisionID; + cursor(0, + collisionID, + generatorID, + getEventInfo(header, Key::pdfParton1Id, int(0)), + getEventInfo(header, Key::pdfParton2Id, int(0)), + getEventInfo(header, Key::pdfCode1, int(0)), + getEventInfo(header, Key::pdfCode2, int(0)), + getEventInfo(header, Key::pdfX1, 0.f), + getEventInfo(header, Key::pdfX2, 0.f), + getEventInfo(header, Key::pdfScale, 1.f), + getEventInfo(header, Key::pdfXF1, 0.f), + getEventInfo(header, Key::pdfXF2, 0.f)); + + return true; +} +//-------------------------------------------------------------------- +bool updateHepMCHeavyIon(const HeavyIonCursor& cursor, + int collisionID, + short generatorID, + o2::dataformats::MCEventHeader const& header, + HepMCUpdate when) +{ + using Key = dataformats::MCInfoKeys; + + if (when == HepMCUpdate::never or + (when != HepMCUpdate::always and // clang + not hasKeys(header, {Key::nCollHard, // format + Key::nPartProjectile, // is + Key::nPartTarget, // so + Key::nColl, // annoying + Key::nCollNNWounded, // . + Key::nCollNWoundedN, // It + Key::nCollNWoundedNwounded, // messes + Key::nSpecProjectileNeutron, // up + Key::nSpecTargetNeutron, // the + Key::nSpecProjectileProton, // clarity + Key::nSpecTargetProton, // of + Key::planeAngle, // the + "eccentricity", // code + Key::sigmaInelNN, // to + Key::centrality}, // noavail + when == HepMCUpdate::anyKey))) { + return false; + } + + int specNeutrons = (getEventInfo(header, Key::nSpecProjectileNeutron, -1) + + getEventInfo(header, Key::nSpecTargetNeutron, -1)); + int specProtons = (getEventInfo(header, Key::nSpecProjectileProton, -1) + + getEventInfo(header, Key::nSpecTargetProton, -1)); + + LOG(verbosity) << "Updating HepMC heavy-ion table w/collisionID " + << collisionID; + cursor(0, + collisionID, + generatorID, + getEventInfo(header, Key::nCollHard, -1), + getEventInfo(header, Key::nPartProjectile, -1), + getEventInfo(header, Key::nPartTarget, -1), + getEventInfo(header, Key::nColl, -1), + getEventInfo(header, Key::nCollNNWounded, -1), + getEventInfo(header, Key::nCollNWoundedN, -1), + getEventInfo(header, Key::nCollNWoundedNwounded, -1), + specNeutrons, + specProtons, + header.GetB(), + getEventInfo(header, Key::planeAngle, header.GetRotZ()), + getEventInfo(header, "eccentricity", 0), + getEventInfo(header, Key::sigmaInelNN, 0.), + getEventInfo(header, Key::centrality, -1)); + + return true; +} +//-------------------------------------------------------------------- +void updateParticle(const ParticleCursor& cursor, + const TrackToIndex& toStore, + int collisionID, + o2::MCTrack const& track, + std::vector const& tracks, + uint8_t flags, + uint32_t weightMask, + uint32_t momentumMask, + uint32_t positionMask) +{ + using o2::mcutils::MCTrackNavigator; + using namespace o2::aod::mcparticle::enums; + using namespace o2::math_utils; + using namespace o2::mcgenstatus; + + auto mapping = [&toStore](int trackNo) { + auto iter = toStore.find(trackNo); + if (iter == toStore.end()) { + return -1; + } + return iter->second; + }; + + auto statusCode = track.getStatusCode().fullEncoding; + auto hepmc = getHepMCStatusCode(track.getStatusCode()); + if (not track.isPrimary()) { + flags |= ProducedByTransport; + statusCode = track.getProcess(); + } + if (MCTrackNavigator::isPhysicalPrimary(track, tracks)) { + flags |= PhysicalPrimary; + } + + int daughters[2] = {-1, -1}; + std::vector mothers; + int id; + if ((id = mapping(track.getMotherTrackId())) >= 0) { + mothers.push_back(id); + } + if ((id = mapping(track.getSecondMotherTrackId())) >= 0) { + mothers.push_back(id); + } + if ((id = mapping(track.getFirstDaughterTrackId())) >= 0) { + daughters[0] = id; + } + if ((id = mapping(track.getFirstDaughterTrackId())) >= 0) { + daughters[1] = id; + } else { + daughters[1] = daughters[0]; + } + if (daughters[0] < 0 and daughters[1] >= 0) { + LOG(error) << "Problematic daughters: " << daughters[0] << " and " + << daughters[1]; + daughters[0] = daughters[1]; + } + if (daughters[0] > daughters[1]) { + std::swap(daughters[0], daughters[1]); + } + + float weight = track.getWeight(); + float pX = float(track.Px()); + float pY = float(track.Py()); + float pZ = float(track.Pz()); + float energy = float(track.GetEnergy()); + float vX = float(track.Vx()); + float vY = float(track.Vy()); + float vZ = float(track.Vz()); + float time = float(track.T()); + + LOG(verbosity) << "McParticle collisionId=" << collisionID << "," + << "status=" << statusCode << "," + << "hepmc=" << hepmc << "," + << "pdg=" << track.GetPdgCode() << "," + << "nMothers=" << mothers.size() << "," + << "daughters=" << daughters[0] << "," + << daughters[1]; + + cursor(0, + collisionID, + track.GetPdgCode(), + statusCode, + flags, + mothers, + daughters, + truncateFloatFraction(weight, weightMask), + truncateFloatFraction(pX, momentumMask), + truncateFloatFraction(pY, momentumMask), + truncateFloatFraction(pZ, momentumMask), + truncateFloatFraction(energy, momentumMask), + truncateFloatFraction(vX, positionMask), + truncateFloatFraction(vY, positionMask), + truncateFloatFraction(vZ, positionMask), + truncateFloatFraction(time, positionMask)); +} +//-------------------------------------------------------------------- +uint32_t updateParticles(const ParticleCursor& cursor, + int collisionID, + std::vector const& tracks, + TrackToIndex& preselect, + uint32_t offset, + bool filter, + bool background, + uint32_t weightMask, + uint32_t momentumMask, + uint32_t positionMask) +{ + using o2::mcutils::MCTrackNavigator; + using namespace o2::aod::mcparticle::enums; + using namespace o2::mcgenstatus; + + // First loop over particles to find out which to store + // TrackToIndex toStore(preselect.begin(), preselect.end()); + // + // Guess we need to modifiy the passed in mapping so that MC labels + // can be set correctly + TrackToIndex& toStore = preselect; + + // Mapping lambda. This maps the track number to the index into + // the table exported. + auto mapping = [&toStore](int trackNo) { + auto iter = toStore.find(trackNo); + if (iter == toStore.end()) { + return -1; + } + return iter->second; + }; + + LOG(verbosity) << "Got a total of " << tracks.size(); + for (int trackNo = tracks.size() - 1; trackNo >= 0; trackNo--) { + auto& track = tracks[trackNo]; + auto hepmc = getHepMCStatusCode(track.getStatusCode()); + if (filter) { + if (toStore.find(trackNo) == toStore.end() and + /* The above test is in-correct. The track may be stored in + the list, but with a negative value. In that case, the + above test will still check mothers and daughters, and + possible store them in the output. This is however how + it is done in current `dev` branch, and so to enable + comparison on closure, we do this test for now. The + correct way it commented out below. */ + // mapping(trackNo) < 0 and + not track.isPrimary() and + not MCTrackNavigator::isPhysicalPrimary(track, tracks) and + not MCTrackNavigator::isKeepPhysics(track, tracks)) { + LOG(verbosity) << "Skipping track " << trackNo << " " << hepmc << " " + << mapping(trackNo) << " " + << track.isPrimary() << " " + << MCTrackNavigator::isPhysicalPrimary(track, tracks) + << " " + << MCTrackNavigator::isKeepPhysics(track, tracks); + continue; + } + } + + // Store this particle. We mark that putting a 1 in the + // `toStore` mapping. This will later on be updated with the + // actual index into the table + toStore[trackNo] = 1; + + // If we're filtering tracks, then also mark mothers and + // daughters(?) to be stored. + if (filter) { + int id; + if ((id = track.getMotherTrackId()) >= 0) { + toStore[id] = 1; + } + if ((id = track.getSecondMotherTrackId()) >= 0) { + toStore[id] = 1; + } + if ((id = track.getFirstDaughterTrackId()) >= 0) { + toStore[id] = 1; + } + if ((id = track.getLastDaughterTrackId()) >= 0) { + toStore[id] = 1; + } + } + } + + // Second loop to set indexes. This is needed to be done before + // we actually update the table, because a particle may point to a + // later particle. + LOG(verbosity) << "Will store " << toStore.size() << " particles"; + size_t index = 0; + + for (size_t trackNo = 0U; trackNo < tracks.size(); trackNo++) { + auto storeIt = mapping(trackNo); + if (storeIt < 0) { + continue; + } + + toStore[trackNo] = offset + index; + index++; + } + + // Make sure we have the right number + assert(index == toStore.size()); + LOG(verbosity) << "Starting index " << offset << ", last index " + << offset + index << " " + << "Selected " << toStore.size() << " tracks out of " + << tracks.size() << " " + << "Collision # " << collisionID; + // Third loop to actually store the particles in the order given + for (size_t trackNo = 0U; trackNo < tracks.size(); trackNo++) { + auto storeIt = mapping(trackNo); + if (storeIt < 0) { + continue; + } + + auto& track = tracks[trackNo]; + auto hepmc = getHepMCStatusCode(track.getStatusCode()); + uint8_t flags = (background ? FromBackgroundEvent : 0); + updateParticle(cursor, + toStore, + collisionID, + track, + tracks, + flags, + weightMask, + momentumMask, + positionMask); + } + LOG(verbosity) << "Return new offset " << offset + index; + return offset + index; +} +} // namespace o2::aodmchelpers + +// Local Variables: +// mode: C++ +// End: diff --git a/Detectors/AOD/src/AODMcProducerWorkflowSpec.cxx b/Detectors/AOD/src/AODMcProducerWorkflowSpec.cxx index 4c1623188e55c..8136bd3e6d268 100644 --- a/Detectors/AOD/src/AODMcProducerWorkflowSpec.cxx +++ b/Detectors/AOD/src/AODMcProducerWorkflowSpec.cxx @@ -13,156 +13,18 @@ #include "AODProducerWorkflow/AODMcProducerWorkflowSpec.h" #include "AODProducerWorkflow/AODProducerHelpers.h" -#include "Framework/AnalysisDataModel.h" #include "Framework/ControlService.h" #include "Framework/DataTypes.h" -#include "Framework/TableBuilder.h" -#include "SimulationDataFormat/MCTrack.h" #include "SimulationDataFormat/MCUtils.h" #include "O2Version.h" #include "TString.h" -#include using namespace o2::framework; -using namespace o2::math_utils::detail; namespace o2::aodmcproducer { -template -void AODMcProducerWorkflowDPL::fillMCParticlesTable(o2::steer::MCKinematicsReader& mcReader, - const MCParticlesCursorType& mcParticlesCursor) -{ - using o2::aodhelpers::Triplet_t; - - int tableIndex = 1; - for (auto& colInfo : mMCColToEvSrc) { // loop over " <-> combined MC col. ID" key pairs - int event = colInfo[2]; - int source = colInfo[1]; - int mcColId = colInfo[0]; - std::vector const& mcParticles = mcReader.getTracks(source, event); - // mark tracks to be stored per event - // loop over stack of MC particles from end to beginning: daughters are stored after mothers - if (mFilterMC) { - for (int particle = mcParticles.size() - 1; particle >= 0; particle--) { - // we store all primary particles == particles put by generator - if (mcParticles[particle].isPrimary() || - o2::mcutils::MCTrackNavigator::isPhysicalPrimary(mcParticles[particle], mcParticles) || - o2::mcutils::MCTrackNavigator::isKeepPhysics(mcParticles[particle], mcParticles)) { - mToStore[Triplet_t(source, event, particle)] = 1; - } else { - continue; - } - // we store mothers and daughters of particles that are to be saved - int mother0 = mcParticles[particle].getMotherTrackId(); - if (mother0 != -1) { - mToStore[Triplet_t(source, event, mother0)] = 1; - } - int mother1 = mcParticles[particle].getSecondMotherTrackId(); - if (mother1 != -1) { - mToStore[Triplet_t(source, event, mother1)] = 1; - } - int daughter0 = mcParticles[particle].getFirstDaughterTrackId(); - if (daughter0 != -1) { - mToStore[Triplet_t(source, event, daughter0)] = 1; - } - int daughterL = mcParticles[particle].getLastDaughterTrackId(); - if (daughterL != -1) { - mToStore[Triplet_t(source, event, daughterL)] = 1; - } - } - // enumerate saved mc particles and their relatives to get mother/daughter relations - for (auto particle = 0U; particle < mcParticles.size(); ++particle) { - auto mapItem = mToStore.find(Triplet_t(source, event, particle)); - if (mapItem != mToStore.end()) { - mapItem->second = tableIndex - 1; - tableIndex++; - } - } - } else { - // if all mc particles are stored, all mc particles will be enumerated - for (auto particle = 0U; particle < mcParticles.size(); ++particle) { - mToStore[Triplet_t(source, event, particle)] = tableIndex - 1; - tableIndex++; - } - } - - // second part: fill survived mc tracks into the AOD table - for (auto particle = 0U; particle < mcParticles.size(); ++particle) { - if (mToStore.find(Triplet_t(source, event, particle)) == mToStore.end()) { - continue; - } - int statusCode = 0; - uint8_t flags = 0; - if (!mcParticles[particle].isPrimary()) { - flags |= o2::aod::mcparticle::enums::ProducedByTransport; // mark as produced by transport - statusCode = mcParticles[particle].getProcess(); - } else { - statusCode = mcParticles[particle].getStatusCode().fullEncoding; - } - if (source == 0) { - flags |= o2::aod::mcparticle::enums::FromBackgroundEvent; // mark as particle from background event - } - if (o2::mcutils::MCTrackNavigator::isPhysicalPrimary(mcParticles[particle], mcParticles)) { - flags |= o2::aod::mcparticle::enums::PhysicalPrimary; // mark as physical primary - } - float weight = mcParticles[particle].getWeight(); - std::vector mothers; - int mcMother0 = mcParticles[particle].getMotherTrackId(); - auto item = mToStore.find(Triplet_t(source, event, mcMother0)); - if (item != mToStore.end()) { - mothers.push_back(item->second); - } - int mcMother1 = mcParticles[particle].getSecondMotherTrackId(); - item = mToStore.find(Triplet_t(source, event, mcMother1)); - if (item != mToStore.end()) { - mothers.push_back(item->second); - } - int daughters[2] = {-1, -1}; // slice - int mcDaughter0 = mcParticles[particle].getFirstDaughterTrackId(); - item = mToStore.find(Triplet_t(source, event, mcDaughter0)); - if (item != mToStore.end()) { - daughters[0] = item->second; - } - int mcDaughterL = mcParticles[particle].getLastDaughterTrackId(); - item = mToStore.find(Triplet_t(source, event, mcDaughterL)); - if (item != mToStore.end()) { - daughters[1] = item->second; - if (daughters[0] < 0) { - LOG(error) << "AOD problematic daughter case observed"; - daughters[0] = daughters[1]; /// Treat the case of first negative label (pruned in the kinematics) - } - } else { - daughters[1] = daughters[0]; - } - if (daughters[0] > daughters[1]) { - std::swap(daughters[0], daughters[1]); - } - auto pX = (float)mcParticles[particle].Px(); - auto pY = (float)mcParticles[particle].Py(); - auto pZ = (float)mcParticles[particle].Pz(); - auto energy = (float)mcParticles[particle].GetEnergy(); - mcParticlesCursor(0, - mcColId, - mcParticles[particle].GetPdgCode(), - statusCode, - flags, - mothers, - daughters, - truncateFloatFraction(weight, mMcParticleW), - truncateFloatFraction(pX, mMcParticleMom), - truncateFloatFraction(pY, mMcParticleMom), - truncateFloatFraction(pZ, mMcParticleMom), - truncateFloatFraction(energy, mMcParticleMom), - truncateFloatFraction((float)mcParticles[particle].Vx(), mMcParticlePos), - truncateFloatFraction((float)mcParticles[particle].Vy(), mMcParticlePos), - truncateFloatFraction((float)mcParticles[particle].Vz(), mMcParticlePos), - truncateFloatFraction((float)mcParticles[particle].T(), mMcParticlePos)); - } - mcReader.releaseTracksForSourceAndEvent(source, event); - } -} - +//------------------------------------------------------------------ void AODMcProducerWorkflowDPL::init(InitContext& ic) { mTimer.Stop(); @@ -172,13 +34,13 @@ void AODMcProducerWorkflowDPL::init(InitContext& ic) mRecoPass = ic.options().get("reco-pass"); mTFNumber = ic.options().get("aod-timeframe-id"); mFilterMC = ic.options().get("filter-mctracks"); - mTruncate = ic.options().get("enable-truncation"); + int truncate = ic.options().get("enable-truncation"); if (mTFNumber == -1L) { LOG(info) << "TFNumber will be obtained from CCDB"; } // set no truncation if selected by user - if (mTruncate != 1) { + if (truncate == 0) { LOG(info) << "Truncation is not used!"; mCollisionPosition = 0xFFFFFFFF; mMcParticleW = 0xFFFFFFFF; @@ -192,112 +54,218 @@ void AODMcProducerWorkflowDPL::init(InitContext& ic) // parse list of sim prefixes into vector mSimPrefix = ic.options().get("mckine-fname"); } + std::string hepmcUpdate = ic.options().get("hepmc-update"); + HepMCUpdate when = (hepmcUpdate == "never" // + ? HepMCUpdate::never // + : hepmcUpdate == "always" // + ? HepMCUpdate::always // + : hepmcUpdate == "all" // + ? HepMCUpdate::allKeys // + : HepMCUpdate::anyKey); + mXSectionUpdate = when; + mPdfInfoUpdate = when; + mHeavyIonUpdate = when; mTimer.Reset(); } +//------------------------------------------------------------------ +void AODMcProducerWorkflowDPL::updateHeader(CollisionCursor& collisionCursor, + XSectionCursor& xSectionCursor, + PdfInfoCursor& pdfInfoCursor, + HeavyIonCursor& heavyIonCursor, + const MCEventHeader& header, + int collisionID, // Index + int bcID, + float time, + short generatorID, + int sourceID) +{ + using aodmchelpers::updateHepMCHeavyIon; + using aodmchelpers::updateHepMCPdfInfo; + using aodmchelpers::updateHepMCXSection; + using aodmchelpers::updateMCCollisions; + + auto genID = updateMCCollisions(collisionCursor, + bcID, + time, + header, + generatorID, + sourceID, + mCollisionPosition); + mXSectionUpdate = (updateHepMCXSection(xSectionCursor, // + collisionID, // + genID, // + header, // + mXSectionUpdate) // + ? HepMCUpdate::always // + : HepMCUpdate::never); + mPdfInfoUpdate = (updateHepMCPdfInfo(pdfInfoCursor, // + collisionID, // + genID, // + header, // + mPdfInfoUpdate) // + ? HepMCUpdate::always // + : HepMCUpdate::never); + mHeavyIonUpdate = (updateHepMCHeavyIon(heavyIonCursor, // + collisionID, // + genID, // + header, // + mHeavyIonUpdate) // + ? HepMCUpdate::always // + : HepMCUpdate::never); +} + +//------------------------------------------------------------------ void AODMcProducerWorkflowDPL::run(ProcessingContext& pc) { mTimer.Start(false); uint64_t tfNumber = mTFNumber; - auto mcCollisionsBuilder = pc.outputs().make(OutputForTable::ref()); - auto mcParticlesBuilder = pc.outputs().make(OutputForTable::ref()); - auto originTableBuilder = pc.outputs().make(OutputForTable::ref()); - - auto mcCollisionsCursor = mcCollisionsBuilder->cursor(); - auto mcParticlesCursor = mcParticlesBuilder->cursor(); - auto originCursor = originTableBuilder->cursor(); - - std::unique_ptr mcReader; - - if (!mEnableEmbed) { - mcReader = std::make_unique(mSimPrefix, steer::MCKinematicsReader::Mode::kMCKine); + using namespace o2::aodmchelpers; + using namespace o2::aodhelpers; + + auto collisionsCursor = createTableCursor(pc); + auto particlesCursor = createTableCursor(pc); + auto originCursor = createTableCursor(pc); + auto xSectionCursor = createTableCursor(pc); + auto pdfInfoCursor = createTableCursor(pc); + auto heavyIonCursor = createTableCursor(pc); + + // --- Create our reader ------------------------------------------- + std::unique_ptr reader; + if (not mEnableEmbed) { + reader = + std::make_unique(mSimPrefix, + MCKinematicsReader::Mode::kMCKine); } else { - mcReader = std::make_unique("collisioncontext.root"); + reader = std::make_unique("collisioncontext.root"); } - // filling mcCollision table + // --- Container of event indexes --------------------------------- + using EventInfo = std::vector>; + EventInfo eventInfo; + + // --- Fill collision and HepMC aux tables ------------------------ // dummy time information - int bcID = 0; float time = 0; - auto updateMCCollisions = [this, bcID, time, &mcCollisionsCursor](dataformats::MCEventHeader const& header, short generatorID, int sourceID) { - bool isValid = false; - int subGeneratorId{-1}; - if (header.hasInfo(o2::mcgenid::GeneratorProperty::SUBGENERATORID)) { - subGeneratorId = header.getInfo(o2::mcgenid::GeneratorProperty::SUBGENERATORID, isValid); - } - float mcColWeight = 1.; - if (header.hasInfo("weight")) { - mcColWeight = header.getInfo("weight", isValid); - } - mcCollisionsCursor(0, - bcID, - o2::mcgenid::getEncodedGenId(header.getInfo(o2::mcgenid::GeneratorProperty::GENERATORID, isValid), sourceID, subGeneratorId), - truncateFloatFraction(header.GetX(), mCollisionPosition), - truncateFloatFraction(header.GetY(), mCollisionPosition), - truncateFloatFraction(header.GetZ(), mCollisionPosition), - truncateFloatFraction(time, mCollisionPosition), - truncateFloatFraction(mcColWeight, mCollisionPosition), - header.GetB()); - }; - - if (!mEnableEmbed) { // simply store all MC events into table + if (not mEnableEmbed) { + // simply store all MC events into table int icol = 0; - int nSources = mcReader->getNSources(); + int nSources = reader->getNSources(); for (int isrc = 0; isrc < nSources; isrc++) { short generatorID = isrc; - int nEvents = mcReader->getNEvents(isrc); + int nEvents = reader->getNEvents(isrc); for (int ievt = 0; ievt < nEvents; ievt++) { - auto& header = mcReader->getMCEventHeader(isrc, ievt); - updateMCCollisions(header, generatorID, isrc); - mMCColToEvSrc.emplace_back(std::vector{icol, isrc, ievt}); + auto& header = reader->getMCEventHeader(isrc, ievt); + updateHeader(collisionsCursor.cursor, + xSectionCursor.cursor, + pdfInfoCursor.cursor, + heavyIonCursor.cursor, + header, + ievt, + ievt, // BC is the same as collision index + time, + generatorID, + isrc); + + eventInfo.emplace_back(std::make_tuple(icol, isrc, ievt)); icol++; } } - } else { // treat embedded events using collisioncontext: injected events will be stored together with background events into the same collisions - int nMCCollisions = mcReader->getDigitizationContext()->getNCollisions(); - const auto& mcRecords = mcReader->getDigitizationContext()->getEventRecords(); - const auto& mcParts = mcReader->getDigitizationContext()->getEventParts(); - for (int icol = 0; icol < nMCCollisions; icol++) { - auto& colParts = mcParts[icol]; + } else { + // treat embedded events using collisioncontext: injected events + // will be stored together with background events into the same + // collisions + int nCollisions = reader->getDigitizationContext()->getNCollisions(); + const auto& records = reader->getDigitizationContext()->getEventRecords(); + const auto& parts = reader->getDigitizationContext()->getEventParts(); + for (int icol = 0; icol < nCollisions; icol++) { + auto& colParts = parts[icol]; auto nParts = colParts.size(); for (auto colPart : colParts) { auto eventID = colPart.entryID; auto sourceID = colPart.sourceID; - // enable embedding: if several colParts exist, then they are saved as one collision + // enable embedding: if several colParts exist, then they are + // saved as one collision if (nParts == 1 || sourceID == 0) { + // Make collision header from first source only short generatorID = sourceID; - auto& header = mcReader->getMCEventHeader(sourceID, eventID); - updateMCCollisions(header, generatorID, sourceID); + auto& header = reader->getMCEventHeader(sourceID, eventID); + + updateHeader(collisionsCursor.cursor, + xSectionCursor.cursor, + pdfInfoCursor.cursor, + heavyIonCursor.cursor, + header, + icol, + icol, // BC is the same as collision index + time, + generatorID, + sourceID); } - mMCColToEvSrc.emplace_back(std::vector{icol, sourceID, eventID}); // point background and injected signal events to one collision + // point background and injected signal events to one collision + eventInfo.emplace_back(std::make_tuple(icol, sourceID, eventID)); } } } - std::sort(mMCColToEvSrc.begin(), mMCColToEvSrc.end(), - [](const std::vector& left, const std::vector& right) { return (left[0] < right[0]); }); - - // filling mc particles table - fillMCParticlesTable(*mcReader, mcParticlesCursor); - - mMCColToEvSrc.clear(); - mToStore.clear(); - - originCursor(0, tfNumber); + // Sort the event information + std::sort(eventInfo.begin(), eventInfo.end(), + [](typename EventInfo::const_reference left, + typename EventInfo::const_reference right) { // + return (std::get<0>(left) < std::get<0>(right)); + }); + + // Loop over available events and update the tracks table + size_t offset = 0; + for (auto& colInfo : eventInfo) { + int event = std::get<2>(colInfo); + int source = std::get<1>(colInfo); + int collisionID = std::get<0>(colInfo); + auto tracks = reader->getTracks(source, event); + + TrackToIndex preselect; + offset = updateParticles(particlesCursor.cursor, + collisionID, + tracks, + preselect, + offset, + mFilterMC, + source == 0, + mMcParticleW, + mMcParticleMom, + mMcParticlePos); + + reader->releaseTracksForSourceAndEvent(source, event); + } - // sending metadata to writer - if (!mIsMDSent) { - TString dataType = "MC"; - TString O2Version = o2::fullVersion(); - TString ROOTVersion = ROOT_RELEASE; - mMetaDataKeys = {"DataType", "Run", "O2Version", "ROOTVersion", "RecoPassName", "AnchorProduction", "AnchorPassName", "LPMProductionTag"}; - mMetaDataVals = {dataType, "3", O2Version, ROOTVersion, mRecoPass, mAnchorProd, mAnchorPass, mLPMProdTag}; - pc.outputs().snapshot(Output{"AMD", "AODMetadataKeys", 0}, mMetaDataKeys); - pc.outputs().snapshot(Output{"AMD", "AODMetadataVals", 0}, mMetaDataVals); + // --- Update the origin and the time-frame ------------------------ + originCursor(tfNumber); + + // --- Sending metadata to writer if not done already -------------- + if (not mIsMDSent) { + TString o2Version = o2::fullVersion(); + std::vector metaDataKeys = {"DataType", + "Run", + "O2Version", + "ROOTVersion", + "RecoPassName", + "AnchorProduction", + "AnchorPassName", + "LPMProductionTag"}; + std::vector metaDataVals = {"MC", + "3", + TString{o2::fullVersion()}, + ROOT_RELEASE, + mRecoPass, + mAnchorProd, + mAnchorPass, + mLPMProdTag}; + pc.outputs().snapshot(Output{"AMD", "AODMetadataKeys", 0}, metaDataKeys); + pc.outputs().snapshot(Output{"AMD", "AODMetadataVals", 0}, metaDataVals); mIsMDSent = true; } @@ -318,10 +286,20 @@ void AODMcProducerWorkflowDPL::endOfStream(EndOfStreamContext&) DataProcessorSpec getAODMcProducerWorkflowSpec() { + using McCollisions = AODMcProducerWorkflowDPL::McCollisions; + using McParticles = AODMcProducerWorkflowDPL::McParticles; + using Origins = AODMcProducerWorkflowDPL::Origins; + using XSections = AODMcProducerWorkflowDPL::XSections; + using PdfInfos = AODMcProducerWorkflowDPL::PdfInfos; + using HeavyIons = AODMcProducerWorkflowDPL::HeavyIons; + std::vector outputs{ - OutputForTable::spec(), - OutputForTable::spec(), - OutputForTable::spec(), + OutputForTable::spec(), + OutputForTable::spec(), + OutputForTable::spec(), + OutputForTable::spec(), + OutputForTable::spec(), + OutputForTable::spec(), OutputSpec{"TFN", "TFNumber"}, OutputSpec{"TFF", "TFFilename"}, OutputSpec{"AMD", "AODMetadataKeys"}, @@ -341,7 +319,11 @@ DataProcessorSpec getAODMcProducerWorkflowSpec() ConfigParamSpec{"reco-pass", VariantType::String, "", {"RecoPassName"}}, ConfigParamSpec{"filter-mctracks", VariantType::Int, 1, {"Store only physical primary MC tracks and their mothers/daughters. 0 -- off, != 0 -- on"}}, ConfigParamSpec{"enable-embedding", VariantType::Int, 0, {"Use collisioncontext.root to process embedded events"}}, - ConfigParamSpec{"mckine-fname", VariantType::String, "o2sim", {"MC kinematics file name prefix: e.g. 'o2sim', 'bkg', 'sgn_1'. Used only if 'enable-embedding' is 0"}}}}; + ConfigParamSpec{"mckine-fname", VariantType::String, "o2sim", {"MC kinematics file name prefix: e.g. 'o2sim', 'bkg', 'sgn_1'. Used only if 'enable-embedding' is 0"}}, + ConfigParamSpec{"hepmc-update", VariantType::String, "always", {"When to update HepMC Aux tables: always - force update, never - never update, all - if all keys are present, any - when any key is present (not valid yet)"}}}}; } } // namespace o2::aodmcproducer +// +// EOF +// diff --git a/Detectors/AOD/src/AODProducerWorkflowSpec.cxx b/Detectors/AOD/src/AODProducerWorkflowSpec.cxx index f1f92dbd01b2e..a236a87370d17 100644 --- a/Detectors/AOD/src/AODProducerWorkflowSpec.cxx +++ b/Detectors/AOD/src/AODProducerWorkflowSpec.cxx @@ -12,6 +12,8 @@ /// @file AODProducerWorkflowSpec.cxx #include "AODProducerWorkflow/AODProducerWorkflowSpec.h" +#include "AODProducerWorkflow/AODMcProducerHelpers.h" +#include "AODProducerWorkflow/AODProducerHelpers.h" #include "DataFormatsEMCAL/TriggerRecord.h" #include "DataFormatsEMCAL/EventHandler.h" #include "DataFormatsFT0/RecPoints.h" @@ -808,6 +810,53 @@ void AODProducerWorkflowDPL::addToFwdTracksTable(FwdTracksCursorType& fwdTracksC } } +//------------------------------------------------------------------ +void AODProducerWorkflowDPL::updateMCHeader(MCCollisionCursor& collisionCursor, + XSectionCursor& xSectionCursor, + PdfInfoCursor& pdfInfoCursor, + HeavyIonCursor& heavyIonCursor, + const MCEventHeader& header, + int collisionID, + int bcID, + float time, + short generatorID, + int sourceID) +{ + using aodmchelpers::updateHepMCHeavyIon; + using aodmchelpers::updateHepMCPdfInfo; + using aodmchelpers::updateHepMCXSection; + using aodmchelpers::updateMCCollisions; + + auto genID = updateMCCollisions(collisionCursor, + bcID, + time, + header, + generatorID, + sourceID, + mCollisionPosition); + mXSectionUpdate = (updateHepMCXSection(xSectionCursor, // + collisionID, // + genID, // + header, // + mXSectionUpdate) + ? HepMCUpdate::always + : HepMCUpdate::never); + mPdfInfoUpdate = (updateHepMCPdfInfo(pdfInfoCursor, // + collisionID, // + genID, // + header, // + mPdfInfoUpdate) + ? HepMCUpdate::always + : HepMCUpdate::never); + mHeavyIonUpdate = (updateHepMCHeavyIon(heavyIonCursor, // + collisionID, // + genID, // + header, + mHeavyIonUpdate) + ? HepMCUpdate::always + : HepMCUpdate::never); +} + void dimensionMCKeepStore(std::vector>>& store, int Nsources, int NEvents) { store.resize(Nsources); @@ -835,9 +884,8 @@ void keepMCParticle(std::vector>>& stor store[source][event][track] = value; } -template void AODProducerWorkflowDPL::fillMCParticlesTable(o2::steer::MCKinematicsReader& mcReader, - MCParticlesCursorType& mcParticlesCursor, + MCParticlesCursor& mcParticlesCursor, const gsl::span& primVer2TRefs, const gsl::span& GIndices, const o2::globaltracking::RecoContainer& data, @@ -918,143 +966,29 @@ void AODProducerWorkflowDPL::fillMCParticlesTable(o2::steer::MCKinematicsReader& keepMCParticle(mToStore, mcTruth.getSourceID(), mcTruth.getEventID(), mcTruth.getTrackID()); } } - int tableIndex = 1; + using namespace aodmchelpers; + using MCTrackNavigator = o2::mcutils::MCTrackNavigator; + + size_t offset = 0; for (auto& colInfo : mcColToEvSrc) { // loop over " <-> combined MC col. ID" key pairs int event = colInfo[2]; int source = colInfo[1]; int mcColId = colInfo[0]; std::vector const& mcParticles = mcReader.getTracks(source, event); - // mark tracks to be stored per event - // loop over stack of MC particles from end to beginning: daughters are stored after mothers - if (mRecoOnly) { - for (int particle = mcParticles.size() - 1; particle >= 0; particle--) { - // we store all primary particles == particles put by generator - if (mcParticles[particle].isPrimary()) { - keepMCParticle(mToStore, source, event, particle); - } else if (o2::mcutils::MCTrackNavigator::isPhysicalPrimary(mcParticles[particle], mcParticles)) { - keepMCParticle(mToStore, source, event, particle); - } else if (o2::mcutils::MCTrackNavigator::isKeepPhysics(mcParticles[particle], mcParticles)) { - keepMCParticle(mToStore, source, event, particle); - } - - // skip treatment if this particle has not been marked during reconstruction - // or based on criteria just above - if (mToStore[source][event].size() > 0 && mToStore[source][event].find(particle) == mToStore[source][event].end()) { - continue; - } - - int mother0 = mcParticles[particle].getMotherTrackId(); - // we store mothers and daughters of particles that are reconstructed - if (mother0 != -1) { - keepMCParticle(mToStore, source, event, mother0); - } - int mother1 = mcParticles[particle].getSecondMotherTrackId(); - if (mother1 != -1) { - keepMCParticle(mToStore, source, event, mother1); - } - int daughter0 = mcParticles[particle].getFirstDaughterTrackId(); - if (daughter0 != -1) { - keepMCParticle(mToStore, source, event, daughter0); - } - int daughterL = mcParticles[particle].getLastDaughterTrackId(); - if (daughterL != -1) { - keepMCParticle(mToStore, source, event, daughterL); - } - } - particleIDsToKeep.clear(); - if (mToStore[source][event].size() > 0) { - LOG(debug) << "The fraction of MC particles kept is " << mToStore[source][event].size() / (1. * mcParticles.size()) << " for source " << source << " and event " << event; - - for (auto& p : mToStore[source][event]) { - particleIDsToKeep.push_back(p.first); - } - std::sort(particleIDsToKeep.begin(), particleIDsToKeep.end()); - for (auto pid : particleIDsToKeep) { - (mToStore[source][event])[pid] = tableIndex - 1; - tableIndex++; - } - } else { - LOG(warn) << "Empty MC event for event id " << event; - } - } else { - // if all mc particles are stored, all mc particles will be enumerated - particleIDsToKeep.clear(); - for (auto particle = 0U; particle < mcParticles.size(); particle++) { - keepMCParticle(mToStore, source, event, particle, tableIndex - 1); - tableIndex++; - particleIDsToKeep.push_back(particle); - } - } + LOG(debug) << "Event=" << event << " source=" << source << " collision=" << mcColId; + auto& preselect = mToStore[source][event]; + + offset = updateParticles(mcParticlesCursor, + mcColId, + mcParticles, + preselect, + offset, + mRecoOnly, + source == 0, // background + mMcParticleW, + mMcParticleMom, + mMcParticlePos); - // second part: fill survived mc tracks into the AOD table - mcParticlesCursor.reserve(particleIDsToKeep.size()); - for (auto particle : particleIDsToKeep) { - int statusCode = 0; - uint8_t flags = 0; - if (!mcParticles[particle].isPrimary()) { - flags |= o2::aod::mcparticle::enums::ProducedByTransport; // mark as produced by transport - statusCode = mcParticles[particle].getProcess(); - } else { - statusCode = mcParticles[particle].getStatusCode().fullEncoding; - } - if (source == 0) { - flags |= o2::aod::mcparticle::enums::FromBackgroundEvent; // mark as particle from background event - } - if (o2::mcutils::MCTrackNavigator::isPhysicalPrimary(mcParticles[particle], mcParticles)) { - flags |= o2::aod::mcparticle::enums::PhysicalPrimary; // mark as physical primary - } - float weight = mcParticles[particle].getWeight(); - std::vector mothers; - int mcMother0 = mcParticles[particle].getMotherTrackId(); - auto item = mToStore[source][event].find(mcMother0); - if (item != mToStore[source][event].end()) { - mothers.push_back(item->second); - } - int mcMother1 = mcParticles[particle].getSecondMotherTrackId(); - item = mToStore[source][event].find(mcMother1); - if (item != mToStore[source][event].end()) { - mothers.push_back(item->second); - } - int daughters[2] = {-1, -1}; // slice - int mcDaughter0 = mcParticles[particle].getFirstDaughterTrackId(); - item = mToStore[source][event].find(mcDaughter0); - if (item != mToStore[source][event].end()) { - daughters[0] = item->second; - } - int mcDaughterL = mcParticles[particle].getLastDaughterTrackId(); - item = mToStore[source][event].find(mcDaughterL); - if (item != mToStore[source][event].end()) { - daughters[1] = item->second; - if (daughters[0] < 0) { - LOG(error) << "AOD problematic daughter case observed"; - daughters[0] = daughters[1]; /// Treat the case of first negative label (pruned in the kinematics) - } - } else { - daughters[1] = daughters[0]; - } - if (daughters[0] > daughters[1]) { - std::swap(daughters[0], daughters[1]); - } - auto pX = (float)mcParticles[particle].Px(); - auto pY = (float)mcParticles[particle].Py(); - auto pZ = (float)mcParticles[particle].Pz(); - auto energy = (float)mcParticles[particle].GetEnergy(); - mcParticlesCursor(mcColId, - mcParticles[particle].GetPdgCode(), - statusCode, - flags, - mothers, - daughters, - truncateFloatFraction(weight, mMcParticleW), - truncateFloatFraction(pX, mMcParticleMom), - truncateFloatFraction(pY, mMcParticleMom), - truncateFloatFraction(pZ, mMcParticleMom), - truncateFloatFraction(energy, mMcParticleMom), - truncateFloatFraction((float)mcParticles[particle].Vx(), mMcParticlePos), - truncateFloatFraction((float)mcParticles[particle].Vy(), mMcParticlePos), - truncateFloatFraction((float)mcParticles[particle].Vz(), mMcParticlePos), - truncateFloatFraction((float)mcParticles[particle].T(), mMcParticlePos)); - } mcReader.releaseTracksForSourceAndEvent(source, event); } } @@ -1775,6 +1709,14 @@ void AODProducerWorkflowDPL::init(InitContext& ic) mZDCTDCMap[ic] = -std::numeric_limits::infinity(); } + std::string hepmcUpdate = ic.options().get("hepmc-update"); + HepMCUpdate when = (hepmcUpdate == "never" ? HepMCUpdate::never : hepmcUpdate == "always" ? HepMCUpdate::always + : hepmcUpdate == "all" ? HepMCUpdate::allKeys + : HepMCUpdate::anyKey); + mXSectionUpdate = when; + mPdfInfoUpdate = when; + mHeavyIonUpdate = when; + mTimer.Reset(); } @@ -1824,6 +1766,8 @@ void AODProducerWorkflowDPL::run(ProcessingContext& pc) LOG(info) << "FOUND " << primVertices.size() << " primary vertices"; + using namespace o2::aodhelpers; + auto bcCursor = createTableCursor(pc); auto cascadesCursor = createTableCursor(pc); auto collisionsCursor = createTableCursor(pc); @@ -1839,6 +1783,9 @@ void AODProducerWorkflowDPL::run(ProcessingContext& pc) auto fwdTrkClsCursor = createTableCursor(pc); auto mcColLabelsCursor = createTableCursor(pc); auto mcCollisionsCursor = createTableCursor(pc); + auto hepmcXSectionsCursor = createTableCursor(pc); + auto hepmcPdfInfosCursor = createTableCursor(pc); + auto hepmcHeavyIonsCursor = createTableCursor(pc); auto mcMFTTrackLabelCursor = createTableCursor(pc); auto mcFwdTrackLabelCursor = createTableCursor(pc); auto mcParticlesCursor = createTableCursor(pc); @@ -1954,6 +1901,8 @@ void AODProducerWorkflowDPL::run(ProcessingContext& pc) std::vector> mcColToEvSrc; if (mUseMC) { + using namespace o2::aodmchelpers; + // filling mcCollision table int nMCCollisions = mcReader->getDigitizationContext()->getNCollisions(); const auto& mcRecords = mcReader->getDigitizationContext()->getEventRecords(); @@ -1974,36 +1923,31 @@ void AODProducerWorkflowDPL::run(ProcessingContext& pc) if (item != bcsMap.end()) { bcID = item->second; } else { - LOG(fatal) << "Error: could not find a corresponding BC ID for MC collision; BC = " << globalBC << ", mc collision = " << iCol; + LOG(fatal) << "Error: could not find a corresponding BC ID " + << "for MC collision; BC = " << globalBC + << ", mc collision = " << iCol; } auto& colParts = mcParts[iCol]; auto nParts = colParts.size(); for (auto colPart : colParts) { auto eventID = colPart.entryID; auto sourceID = colPart.sourceID; - // enable embedding: if several colParts exist, then they are saved as one collision + // enable embedding: if several colParts exist, then they are + // saved as one collision if (nParts == 1 || sourceID == 0) { // FIXME: // use generators' names for generatorIDs (?) auto& header = mcReader->getMCEventHeader(sourceID, eventID); - bool isValid = false; - int subGeneratorId{-1}; - if (header.hasInfo(o2::mcgenid::GeneratorProperty::SUBGENERATORID)) { - subGeneratorId = header.getInfo(o2::mcgenid::GeneratorProperty::SUBGENERATORID, isValid); - } - isValid = false; - float mcColWeight = 1.; - if (header.hasInfo("weight")) { - mcColWeight = header.getInfo("weight", isValid); - } - mcCollisionsCursor(bcID, - o2::mcgenid::getEncodedGenId(header.getInfo(o2::mcgenid::GeneratorProperty::GENERATORID, isValid), sourceID, subGeneratorId), - truncateFloatFraction(header.GetX(), mCollisionPosition), - truncateFloatFraction(header.GetY(), mCollisionPosition), - truncateFloatFraction(header.GetZ(), mCollisionPosition), - truncateFloatFraction(time, mCollisionPosition), - truncateFloatFraction(mcColWeight, mCollisionPosition), - header.GetB()); + updateMCHeader(mcCollisionsCursor.cursor, + hepmcXSectionsCursor.cursor, + hepmcPdfInfosCursor.cursor, + hepmcHeavyIonsCursor.cursor, + header, + iCol, + bcID, + time, + 0, + sourceID); } mcColToEvSrc.emplace_back(std::vector{iCol, sourceID, eventID}); // point background and injected signal events to one collision } @@ -2273,7 +2217,7 @@ void AODProducerWorkflowDPL::run(ProcessingContext& pc) timer.Start(); // filling mc particles table fillMCParticlesTable(*mcReader, - mcParticlesCursor, + mcParticlesCursor.cursor, primVer2TRefs, primVerGIs, recoData, @@ -2577,7 +2521,7 @@ bool AODProducerWorkflowDPL::propagateTrackToPV(o2::track::TrackParametrizationW dcaInfo.set(999.f, 999.f, 999.f, 999.f, 999.f); o2::dataformats::VertexBase v = mVtx.getMeanVertex(colID < 0 ? 0.f : data.getPrimaryVertex(colID).getZ()); return o2::base::Propagator::Instance()->propagateToDCABxByBz(v, trackPar, 2.f, mMatCorr, &dcaInfo); -}; +} void AODProducerWorkflowDPL::extrapolateToCalorimeters(TrackExtraInfo& extraInfoHolder, const o2::track::TrackPar& track) { @@ -2596,7 +2540,8 @@ void AODProducerWorkflowDPL::extrapolateToCalorimeters(TrackExtraInfo& extraInfo SPHOS, // 12 PHOS only SPHOS | SEMCAL, SPHOS | SEMCAL, SPHOS | SEMCAL, // 13:15 PHOS & DCAL SEMCAL, // 16 DCAL only - SNONE}; // 17 + SNONE // 17 + }; o2::track::TrackPar outTr{track}; auto prop = o2::base::Propagator::Instance(); @@ -2964,6 +2909,9 @@ DataProcessorSpec getAODProducerWorkflowSpec(GID::mask_t src, bool enableSV, boo OutputForTable::spec(), OutputForTable::spec(), OutputForTable::spec(), + OutputForTable::spec(), + OutputForTable::spec(), + OutputForTable::spec(), OutputForTable::spec(), OutputForTable::spec(), OutputForTable::spec(), @@ -3017,6 +2965,7 @@ DataProcessorSpec getAODProducerWorkflowSpec(GID::mask_t src, bool enableSV, boo ConfigParamSpec{"ctpreadout-create", VariantType::Int, 0, {"Create CTP digits from detector readout and CTP inputs. !=1 -- off, 1 -- on"}}, ConfigParamSpec{"emc-select-leading", VariantType::Bool, false, {"Flag to select if only the leading contributing particle for an EMCal cell should be stored"}}, ConfigParamSpec{"propagate-tracks", VariantType::Bool, false, {"Propagate tracks (not used for secondary vertices) to IP"}}, + ConfigParamSpec{"hepmc-update", VariantType::String, "always", {"When to update HepMC Aux tables: always - force update, never - never update, all - if all keys are present, any - when any key is present (not valid yet)"}}, ConfigParamSpec{"propagate-muons", VariantType::Bool, false, {"Propagate muons to IP"}}, ConfigParamSpec{"trackqc-fraction", VariantType::Float, float(0.1), {"Fraction of tracks to QC"}}, ConfigParamSpec{"trackqc-NTrCut", VariantType::Int64, 4L, {"Minimal length of the track - in amount of tracklets"}}, diff --git a/Detectors/AOD/src/aod-mc-producer-workflow.cxx b/Detectors/AOD/src/aod-mc-producer-workflow.cxx index 0bfd1ff021272..4f32bacac6fef 100644 --- a/Detectors/AOD/src/aod-mc-producer-workflow.cxx +++ b/Detectors/AOD/src/aod-mc-producer-workflow.cxx @@ -41,3 +41,6 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) wf.emplace_back(o2::aodmcproducer::getAODMcProducerWorkflowSpec()); return std::move(wf); } +// +// EOF +// diff --git a/Generators/include/Generators/AODToHepMC.h b/Generators/include/Generators/AODToHepMC.h index c3e16f694f746..9427c2336d213 100644 --- a/Generators/include/Generators/AODToHepMC.h +++ b/Generators/include/Generators/AODToHepMC.h @@ -24,6 +24,7 @@ #include #include #include +#define AODTOHEPMC_WITH_HEAVYION namespace o2 { @@ -390,6 +391,11 @@ struct AODToHepMC { * configurables. */ virtual void init(); + /** + * Call before starting to process an event. This clears the + * current event and internal data structures. + */ + virtual void startEvent(); /** * Process the collision header and tracks * @@ -407,8 +413,17 @@ struct AODToHepMC { */ virtual void process(Header const& collision, XSections const& xsections, - PdfInfos const& pdfs, - HeavyIons const& heavyions); + PdfInfos const& pdfs +#ifdef AODTOHEPMC_WITH_HEAVYION + , + HeavyIons const& heavyions +#endif + ); + /** + * Call after process an. Thisf finalises the event and optionally + * outputs to dump. + */ + virtual void endEvent(); /** * End of run - closes output file if enabled. This is called via * specialisation of o2::framework::OutputManager. diff --git a/Generators/src/AODToHepMC.cxx b/Generators/src/AODToHepMC.cxx index 8d508e24ee1cb..4524bf7f2a60e 100644 --- a/Generators/src/AODToHepMC.cxx +++ b/Generators/src/AODToHepMC.cxx @@ -35,15 +35,20 @@ void AODToHepMC::init() << " Output precision: " << mPrecision; } // ------------------------------------------------------------------- -void AODToHepMC::process(Header const& collision, - Tracks const& tracks) +void AODToHepMC::startEvent() { - LOG(debug) << "--- Processing track information"; + LOG(debug) << ">>> Starting an event"; mBeams.clear(); mOrphans.clear(); mParticles.clear(); mVertices.clear(); mEvent.clear(); +} +// ------------------------------------------------------------------- +void AODToHepMC::process(Header const& collision, + Tracks const& tracks) +{ + LOG(debug) << "--- Processing track information"; makeHeader(collision); makeParticles(tracks); @@ -52,13 +57,29 @@ void AODToHepMC::process(Header const& collision, // ------------------------------------------------------------------- void AODToHepMC::process(Header const& collision, XSections const& xsections, - PdfInfos const& pdfs, - HeavyIons const& heavyions) + PdfInfos const& pdfs +#ifdef AODTOHEPMC_WITH_HEAVYION + , + HeavyIons const& heavyions +#endif +) { LOG(debug) << "--- Processing auxiliary information"; makeXSection(xsections); makePdfInfo(pdfs); +#ifdef AODTOHEPMC_WITH_HEAVYION makeHeavyIon(heavyions, collision); +#endif +} +// ------------------------------------------------------------------- +void AODToHepMC::endEvent() +{ + LOG(debug) << "<<< an event"; + if (not mWriter) { + return; + } + // If we have a writer, then dump event to output file + mWriter->write_event(mEvent); } // =================================================================== void AODToHepMC::makeEvent(Header const& collision, @@ -105,10 +126,6 @@ void AODToHepMC::makeEvent(Header const& collision, } // 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) @@ -134,6 +151,7 @@ void AODToHepMC::makeHeader(Header const& header) // ------------------------------------------------------------------- void AODToHepMC::makeXSection(XSections const& xsections) { + LOG(debug) << "--- Process cross-section information"; if (not mCrossSec) { // If we do not have a cross-sections object, create it mCrossSec = std::make_shared(); @@ -144,6 +162,7 @@ void AODToHepMC::makeXSection(XSections const& xsections) if (xsections.size() <= 0) { // If we have no info, skip the rest + LOG(warning) << "??? No input cross-section"; return; } @@ -156,6 +175,7 @@ void AODToHepMC::makeXSection(XSections const& xsections) // ------------------------------------------------------------------- void AODToHepMC::makePdfInfo(PdfInfos const& pdfs) { + LOG(debug) << "--- Process PDF information"; if (not mPdf) { // If we do not have a Parton Distribution Function object, create it mPdf = std::make_shared(); @@ -166,6 +186,7 @@ void AODToHepMC::makePdfInfo(PdfInfos const& pdfs) if (pdfs.size() <= 0) { // If we have no PDF info, skip the rest + LOG(warning) << "??? No input pdf-info, skipping"; return; } @@ -184,6 +205,7 @@ void AODToHepMC::makePdfInfo(PdfInfos const& pdfs) void AODToHepMC::makeHeavyIon(HeavyIons const& heavyions, Header const& header) { + LOG(debug) << "--- Process input heavy-ion header"; if (not mIon) { // Generate heavy ion element if it doesn't exist mIon = std::make_shared(); @@ -215,6 +237,7 @@ void AODToHepMC::makeHeavyIon(HeavyIons const& heavyions, if (heavyions.size() <= 0) { // If we have no heavy-ion information, skip the rest + LOG(warning) << "??? No input heavy-ion header, skipping"; return; } diff --git a/run/CMakeLists.txt b/run/CMakeLists.txt index c8efa3ea56f92..58e5d9e5af402 100644 --- a/run/CMakeLists.txt +++ b/run/CMakeLists.txt @@ -125,8 +125,11 @@ o2_add_executable(hepmc-publisher o2_add_executable(mctracks-to-aod COMPONENT_NAME sim - SOURCES o2sim_mctracks_to_aod.cxx + SOURCES o2sim_mctracks_to_aod.cxx ../Detectors/AOD/src/AODMcProducerHelpers.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::SimulationDataFormat) +target_include_directories(O2exe-sim-mctracks-to-aod + PRIVATE + ../Detectors/AOD/include) o2_add_executable(mc-to-hepmc COMPONENT_NAME aod diff --git a/run/o2aod_mc_to_hepmc.cxx b/run/o2aod_mc_to_hepmc.cxx index 091826d8335a7..3b097fb484d42 100644 --- a/run/o2aod_mc_to_hepmc.cxx +++ b/run/o2aod_mc_to_hepmc.cxx @@ -11,10 +11,34 @@ /** @author Christian Holm Christensen */ -#include #include #include #include +// #define HEPMC_PROCESS_AUX + +#ifndef HEPMC_PROCESS_AUX +//-------------------------------------------------------------------- +using o2::framework::ConfigParamKind; +using o2::framework::ConfigParamSpec; + +// ------------------------------------------------------------------- +void customize(std::vector& workflowOptions) +{ + using o2::framework::VariantType; + + workflowOptions.emplace_back(ConfigParamSpec{"hepmc-aux", // + VariantType::Bool, // + false, // + {"Also process auxiliary " // + "HepMC tables"}, + ConfigParamKind::kProcessFlag}); +} +#endif + +//-------------------------------------------------------------------- +// This _must_ be included after our "customize" function above, or +// that function will not be taken into account. +#include //-------------------------------------------------------------------- /** Task to convert AOD MC tables into HepMC event structure @@ -74,18 +98,84 @@ struct Task1 { void process(Header const& collision, XSections const& xsections, PdfInfos const& pdfs, +#ifdef AODTOHEPMC_WITH_HEAVYION HeavyIons const& heavyions, +#endif Tracks const& tracks) { LOG(debug) << "=== Processing everything ==="; + mConverter.startEvent(); mConverter.process(collision, xsections, - pdfs, - heavyions); + pdfs +#ifdef AODTOHEPMC_WITH_HEAVYION + , + heavyions +#endif + ); mConverter.process(collision, tracks); + mConverter.endEvent(); } }; +//-------------------------------------------------------------------- +/** + * Same as Task1 above, except only header and tracks are processed. + * + * - @c o2::aod::McCollisions + * - @c o2::aod::McParticles + * + */ +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(); + } + /** 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 ==="; + mConverter.startEvent(); + mConverter.process(collision, tracks); + mConverter.endEvent(); + } +}; + +#ifdef HEPMC_PROCESS_AUX //-------------------------------------------------------------------- /** * Ideally, this application should work with the case where only @@ -113,6 +203,12 @@ struct Task1 { * - : TFN/TFNumber/0 * @endverbatim * + * Or + * + * @verbatim + * InputRecord::get: no input with binding HepMCHeavyIons found. Available inputs: McCollisions, McParticles + * @endverbatim + * * Interstingly, the application @c o2-sim-mcevent-to-aod works fine * on its own, e.g., like * @@ -125,8 +221,20 @@ struct Task1 { * @endverbatim * * works fine. + * + * Actually, it is not likely that this will ever work. The various + * process are done out of sync. That is, first all input events of + * the timeframe are passed to the regular `process` method - i.e., + * tracks and collision headers are processed. Then all input events + * of the timeframe are passed to the optional `processAux` method - + * i.e., auxiliary tables and collision headers. + * + * That means that we cannot correlate the tracks and aux tables into + * one event, which is what we need to format a proper HepMC + * event. The reason why it could work in the above example is because + * we only process one timeframe at a time. */ -struct Task2 { +struct Task3 { /** Alias the converter type */ using Converter = o2::eventgen::AODToHepMC; @@ -207,15 +315,16 @@ struct Task2 { * command line argument (@c --hepmc-aux) rather than to rely on an * auto-generated name (would be @ --processAux). */ - decltype(o2::framework::ProcessConfigurable{&Task2::processAux, + decltype(o2::framework::ProcessConfigurable{&Task3::processAux, "hepmc-aux", false, "Process auxilirary " "information"}) - doAux = o2::framework::ProcessConfigurable{&Task2::processAux, + doAux = o2::framework::ProcessConfigurable{&Task3::processAux, "hepmc-aux", false, "Process auxilirary " "information"}; }; +#endif //-------------------------------------------------------------------- using WorkflowSpec = o2::framework::WorkflowSpec; @@ -229,9 +338,21 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfg) using o2::framework::adaptAnalysisTask; // Task1: One entry: header, tracks, auxiliary - // Task2: Two entry: header, tracks, and auxiliary + // Task2: One entry: header, tracks + // Task3: Two entry: header, tracks, and auxiliary +#ifndef HEPMC_PROCESS_AUX + if (cfg.options().get("hepmc-aux")) { + LOG(info) << "Creating full o2-aod-to-hepmc processor"; + return WorkflowSpec{ + adaptAnalysisTask(cfg, TaskName{"o2-aod-to-hepmc"})}; + } + LOG(info) << "Creating minimal o2-aod-to-hepmc processor"; + return WorkflowSpec{ + adaptAnalysisTask(cfg, TaskName{"o2-aod-to-hepmc"})}; +#else return WorkflowSpec{ - adaptAnalysisTask(cfg, TaskName{"o2-aod-mc-to-hepmc"})}; + adaptAnalysisTask(cfg, TaskName{"o2-aod-mc-to-hepmc"})}; +#endif } // // EOF diff --git a/run/o2sim_mctracks_to_aod.cxx b/run/o2sim_mctracks_to_aod.cxx index 355eade881f79..c63063a721d43 100644 --- a/run/o2sim_mctracks_to_aod.cxx +++ b/run/o2sim_mctracks_to_aod.cxx @@ -9,182 +9,149 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -#include "Framework/AnalysisTask.h" -#include "Framework/runDataProcessing.h" -#include "SimulationDataFormat/MCTrack.h" -#include "SimulationDataFormat/MCEventHeader.h" -#include "SimulationDataFormat/MCUtils.h" -#include "Framework/AnalysisDataModel.h" -#include "SimulationDataFormat/InteractionSampler.h" +// #include "O2RivetExporter.h" +#include "../Detectors/AOD/include/AODProducerWorkflow/AODMcProducerHelpers.h" +#include +#include +#include -using namespace o2::framework; +template +using Configurable = o2::framework::Configurable; -struct MctracksToAod { - Produces mcollisions; - Produces mcparticles; +struct Task { + /** @{ + @name Types used */ + using Collisions = o2::aod::McCollisions; + // using Particles = o2::aod::McParticles; // + using Particles = o2::aod::StoredMcParticles_001; + using XSections = o2::aod::HepMCXSections; + using PdfInfos = o2::aod::HepMCPdfInfos; + using HeavyIons = o2::aod::HepMCHeavyIons; + using InteractionSampler = o2::steer::InteractionSampler; + /** @} */ - // optional extra info from HepMC (if found in header) - Produces hepmcxsections; - Produces hepmcpdfinfos; - Produces hepmcheavyions; + /** @{ + @name Produced data */ + /** Collision header */ + o2::framework::Produces mCollisions; + /** Particles in collision */ + o2::framework::Produces mParticles; + /** Cross-section information */ + o2::framework::Produces mXSections; + /** Parton-distribution-function information */ + o2::framework::Produces mPdfInfos; + /** Heavy-ion colllision information */ + o2::framework::Produces mHeavyIons; + /** @} */ + /** @{ + @name Configurable parameters */ + Configurable IR{"interaction-rate", 100.f, + "Interaction rate to simulate"}; + Configurable filt{"filter-mctracks", false, + "Filter tracks"}; + /** @} */ - Configurable IR{"interaction-rate", 100.f, "Interaction rate to simulate"}; - - uint64_t timeframe = 0; - o2::steer::InteractionSampler sampler; + /** Number of timeframes */ + uint64_t mTimeFrame = 0; + /** Interaction simulation */ + InteractionSampler mSampler; + /** Whether to filter tracks */ + bool mFilter; + /** Initialize */ void init(o2::framework::InitContext& /*ic*/) { - sampler.setInteractionRate(IR); - sampler.setFirstIR({0, 0}); - sampler.init(); + mSampler.setInteractionRate(IR); + mSampler.setFirstIR({0, 0}); + mSampler.init(); + mFilter = filt; } + /** Run the conversion */ void run(o2::framework::ProcessingContext& pc) { - auto Nparts = pc.inputs().getNofParts(0); - auto Nparts_verify = pc.inputs().getNofParts(1); - if (Nparts != Nparts_verify) { - LOG(warn) << "Mismatch between number of MC headers and number of track vectors: " << Nparts << " != " << Nparts_verify << ", shipping the empty timeframe"; + LOG(info) << "=== Running extended MC AOD exporter ==="; + using namespace o2::aodmchelpers; + using McHeader = o2::dataformats::MCEventHeader; + using McTrack = o2::MCTrack; + using McTracks = std::vector; + + auto nParts = pc.inputs().getNofParts(0); + auto nPartsVerify = pc.inputs().getNofParts(1); + if (nParts != nPartsVerify) { + LOG(warn) << "Mismatch between number of MC headers and " + << "number of track vectors: " << nParts + << " != " << nPartsVerify + << ", shipping the empty timeframe"; return; } // TODO: include BC simulation auto bcCounter = 0UL; - int trackCounter = 0; - for (auto i = 0U; i < Nparts; ++i) { - auto record = sampler.generateCollisionTime(); - auto mcheader = pc.inputs().get("mcheader", i); - auto mctracks = pc.inputs().get>("mctracks", i); - - mcollisions(bcCounter++, // bcId - 0, // generatorId - mcheader->GetX(), - mcheader->GetY(), - mcheader->GetZ(), - record.timeInBCNS * 1.e-3, // ns to ms - 1., // weight - mcheader->GetB()); - bool valid; - if (mcheader->hasInfo("Accepted")) { - auto ptHard = mcheader->hasInfo("PtHard") ? mcheader->getInfo("PtHard", valid) : 0.f; - auto nMPI = mcheader->hasInfo("MPI") ? mcheader->getInfo("MPI", valid) : -1; - auto processId = mcheader->hasInfo("ProcessId") ? mcheader->getInfo("ProcessId", valid) : -1; - hepmcxsections( - mcollisions.lastIndex(), - 0, - mcheader->getInfo("Accepted", valid), - mcheader->getInfo("Attempted", valid), - mcheader->getInfo("XsectGen", valid), - mcheader->getInfo("XsectErr", valid), - ptHard, - nMPI, - processId); - } - if (mcheader->hasInfo("Id1")) { - hepmcpdfinfos( - mcollisions.lastIndex(), - 0, - mcheader->getInfo("Id1", valid), - mcheader->getInfo("Id2", valid), - mcheader->getInfo("PdfId1", valid), - mcheader->getInfo("PdfId2", valid), - mcheader->getInfo("X1", valid), - mcheader->getInfo("X2", valid), - mcheader->getInfo("scale", valid), - mcheader->getInfo("Pdf1", valid), - mcheader->getInfo("Pdf2", valid)); - } - if (mcheader->hasInfo("NcollHard")) { - hepmcheavyions( - mcollisions.lastIndex(), - 0, - mcheader->getInfo("NcollHard", valid), - mcheader->getInfo("NpartProj", valid), - mcheader->getInfo("NpartTarg", valid), - mcheader->getInfo("Ncoll", valid), - mcheader->getInfo("NNwoundedCollisions", valid), - mcheader->getInfo("NwoundedNCollisions", valid), - mcheader->getInfo("NwoundedNwoundedCollisions", valid), - mcheader->getInfo("SpectatorNeutrons", valid), - mcheader->getInfo("SpectatorProtons", valid), - mcheader->getInfo("ImpactParameter", valid), - mcheader->getInfo("EventPlaneAngle", valid), - mcheader->getInfo("Eccentricity", valid), - mcheader->getInfo("SigmaInelNN", valid), - mcheader->getInfo("Centrality", valid)); - } - - for (auto& mctrack : mctracks) { - std::vector mothers; - int daughters[2]; - - if (mctrack.getMotherTrackId() >= 0) { - mothers.push_back(mctrack.getMotherTrackId() + trackCounter); - } - if (mctrack.getSecondMotherTrackId() >= 0) { - mothers.push_back(mctrack.getSecondMotherTrackId() + trackCounter); - } - daughters[0] = -1; - daughters[1] = -1; - if (mctrack.getFirstDaughterTrackId() >= 0 && mctrack.getLastDaughterTrackId() >= 0) { - daughters[0] = mctrack.getFirstDaughterTrackId() + trackCounter; - daughters[1] = mctrack.getLastDaughterTrackId() + trackCounter; - } else if (mctrack.getFirstDaughterTrackId() >= 0) { - daughters[0] = mctrack.getFirstDaughterTrackId() + trackCounter; - daughters[1] = mctrack.getLastDaughterTrackId() + trackCounter; - } - int PdgCode = mctrack.GetPdgCode(); - int statusCode = 0; - float weight = mctrack.getWeight(); - float px = mctrack.Px(); - float py = mctrack.Py(); - float pz = mctrack.Pz(); - float e = mctrack.GetEnergy(); - float x = mctrack.GetStartVertexCoordinatesX(); - float y = mctrack.GetStartVertexCoordinatesY(); - float z = mctrack.GetStartVertexCoordinatesZ(); - float t = mctrack.GetStartVertexCoordinatesT(); - int flags = 0; - if (!mctrack.isPrimary()) { - flags |= o2::aod::mcparticle::enums::ProducedByTransport; // mark as produced by transport - statusCode = mctrack.getProcess(); - } else { - statusCode = mctrack.getStatusCode().fullEncoding; - } - if (o2::mcutils::MCTrackNavigator::isPhysicalPrimary(mctrack, mctracks)) { - flags |= o2::aod::mcparticle::enums::PhysicalPrimary; // mark as physical primary - } - mcparticles(mcollisions.lastIndex(), // collisionId, - PdgCode, - statusCode, - flags, - mothers, - daughters, - weight, - px, - py, - pz, - e, - x, - y, - z, - t); - } - trackCounter = trackCounter + mctracks.size(); + size_t offset = 0; + for (auto i = 0U; i < nParts; ++i) { + LOG(info) << "--- Loop over " << nParts << " parts ---"; + + auto record = mSampler.generateCollisionTime(); + auto header = pc.inputs().get("mcheader", i); + auto tracks = pc.inputs().get("mctracks", i); + + LOG(info) << "Updating collision table"; + auto genID = updateMCCollisions(mCollisions.cursor, + bcCounter, + record.timeInBCNS * 1.e-3, + *header, + 0, + i); + + LOG(info) << "Updating HepMC tables"; + updateHepMCXSection(mXSections.cursor, bcCounter, genID, *header); + updateHepMCPdfInfo(mPdfInfos.cursor, bcCounter, genID, *header); + updateHepMCHeavyIon(mHeavyIons.cursor, bcCounter, genID, *header); + + LOG(info) << "Updating particles table"; + TrackToIndex preselect; + offset = updateParticles(mParticles.cursor, + bcCounter, + tracks, + preselect, + offset, + mFilter, + false); + + LOG(info) << "Increment BC counter"; + bcCounter++; } - ++timeframe; + using o2::framework::Lifetime; + using o2::framework::Output; + + ++mTimeFrame; pc.outputs().snapshot(Output{"TFF", "TFFilename", 0}, ""); - pc.outputs().snapshot(Output{"TFN", "TFNumber", 0}, timeframe); + pc.outputs().snapshot(Output{"TFN", "TFNumber", 0}, mTimeFrame); } }; +using WorkflowSpec = o2::framework::WorkflowSpec; +using TaskName = o2::framework::TaskName; +using DataProcessorSpec = o2::framework::DataProcessorSpec; +using ConfigContext = o2::framework::ConfigContext; +using InputSpec = o2::framework::InputSpec; +using OutputSpec = o2::framework::OutputSpec; +using Lifetime = o2::framework::Lifetime; + WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { - auto dSpec = adaptAnalysisTask(cfgc); - dSpec.inputs.emplace_back("mctracks", "MC", "MCTRACKS", 0., Lifetime::Timeframe); - dSpec.inputs.emplace_back("mcheader", "MC", "MCHEADER", 0., Lifetime::Timeframe); - dSpec.outputs.emplace_back("TFF", "TFFilename"); - dSpec.outputs.emplace_back("TFN", "TFNumber"); + using o2::framework::adaptAnalysisTask; - return {dSpec}; + auto spec = adaptAnalysisTask(cfgc); + spec.inputs.emplace_back("mctracks", "MC", "MCTRACKS", 0., + Lifetime::Timeframe); + spec.inputs.emplace_back("mcheader", "MC", "MCHEADER", 0., + Lifetime::Timeframe); + spec.outputs.emplace_back("TFF", "TFFilename"); + spec.outputs.emplace_back("TFN", "TFNumber"); + + return {spec}; } +// +// EOF +//