From f87bad43db186765c24b22f1b53b6bbb51c30225 Mon Sep 17 00:00:00 2001 From: Christian Holm Christensen Date: Fri, 19 Jan 2024 10:17:31 +0100 Subject: [PATCH 1/2] Refactor MC part of AOD producers The MC part of the AOD producer workflow `o2::aodproducer::AODProducerWorkflowDPL` and `o2::aodmcproducer::AODMcProducerWorkflowDPL` is refactored to use functions from namespace `o2::aodmchelpers`. The helpers are - `updateMCCollisions` which takes in the `MCEventHeader` and writes to the `o2::aod::McCollisions` table. - `updateHepMCXSection` which takes the `MCEventHeader` and writes to the `o2::aodHepMCSections` table. This uses the predefined key constants as defined in `o2::dataformats::MCInfoKeys` - `updateHepMCPdfInfo` similar to `updateHepMCXSection` above - `updateHepMCHeavyIon` similar to `updateHepMCXSection` above - `updateParticle` uses information from an `o2::MCTrack` and writes it to the `o2::aod::McParticles` table - `updateParticles` loops over `o2::MCTrack` objects and calls `updateParticle`. These functions, in particular `updateHepMC...` uses the functions - `hasKeys` which checks if the `MCEventHeader` has any or all of the keys queried. - `getEventInfo` gets auxiliary information from the `MCEventHeader` or a default value. For the `o2::aod::HepMC...` tables: Depending on the policy parameter passed to the `updateHepMC...` functons, these tables may or may not be updated. - If the policy is `HepMCUpdate::never` then the tables are never updated. - If the policy is `HepMCUpdate::always` then the tables are _always_ updated, possibly with default values. - If the policy is `HepMCUpdate::anyKey` (default) or `HepMCUpdate::allKeys`, then the decision of what to do is taken on the first event seen. - If the policy is `HepMCUpdate::anyKey`, then if _any_ of the needed keys are present, then updating will be enabled for this and _all_ subsequent events. - If the policy is `HepMCUpdate::allKeys`, then if _all_ of the needed keys are present, then updating will be enabled for this and _all_ subsequent events. Note that the availability of keys is _not_ checked after the first event. That means, if the criteria isn't met on the first event, then the tables will _never_ be update (as if the policy was `HepMCUpdate::never`). On the other hand, if the criteria was met, than the tables _will_ be update an all events (as if the policy was `HepMCUpdate::always`). Note the slightly tricky template `TableCursor` which allows us to define a type that correponds to a table curser (which is really a lambda). This template could be moved to `AnalysisDataFormats.h` or the like. The applications `o2-aod-producer-workflow` and `o2-aod-mc-producer-workflow` have been updated (via their respective implementation classes) to use these tools, thus unifying how the MC information is propagated to AODs. The utility `o2-sim-mctracks-to-aod` (`run/o2sim_mctracks_to_aod.cxx`) has _also_ been updated to use these common tools. Both `o2-aod-mc-producer-workflow` and `o2-sim-mctracks-to-aod` has been tested extensively. `o2-aod-producer-workflow` has _not_ been tested since it is not clear to me how to set-up such a test with limited resources. However, since the changes _only_ effect the MC part, and that part is now common between the two `o2-aod-{,mc-}producer-workflow` applications, I believe there is no reason to think that it wouldn't work. **Some important (late) bug fixes** Since the commits are squashed into one, I give some more details here on some later commits. For future reference. **HepMC aux tables** I found two problems with the HepMC aux tables after looking at the data a little deeper - I used the BC identifier as the collision index field in the HepMC aux tables. This happened because I took my clues from the MC producer. The MC producer does not generate a BC ID - even though it sort of looked like it - the BC ID in the MC producer is essentially an index. I've fixed this by passing the collision index (called `collisionID` most places) to the relevant member functions and the table updates. - The producers were set up so that HepMC aux tables would _only_ be written if the input had the corresponding data. If a model gave the Cross-section, but not PDF nor Heavy-Ion information, then only the cross-section table would be populated. Pythia, for example, gives the cross-section and PDF information, but no Heavy-Ion information. All three tables would be produced, but some may not have any entries. Later on, when we want to process the events (by Rivet f.ex.), we like to access as much HepMC aux information as possible so we may build the most complete HepMC event possible. Thus, we would like to read in - MC Collisions - MC particles - 3 HepMC aux However, if one or more of the AOD trees of the 3 HepMC aux tables had no entries, it will cause the producer to crash ``` libO2FrameworkAnalysisSupport.so: std::function::operator()(o2::framework::TreeToTable&) const libO2FrameworkAnalysisSupport.so: o2::framework::LifetimeHolder::release() libO2FrameworkAnalysisSupport.so: o2::framework::LifetimeHolder::~LifetimeHolder() libO2FrameworkAnalysisSupport.so: o2::framework::DataInputDescriptor::readTree(o2::framework::DataAllocator&, o2::header::DataHeader, int, int, std::__cxx11::basic_string, std::allocator >, unsigned long&, unsigned long&) ``` I cannot quite figure out why that happens, but I guess it's a problem triggered by `TreeToTable` or the call back on that object. The (temporary) solution to the above problem is to set the update policy on the HepMC aux tables to be `always`. That means we will always write the tables and give them entries. The real solution to the problem will be to fix `TreeToTable` or what ever is causing the above `SIGSEGV`. **MC track lables** To get MC labels correct, the member `AODProducerWorkflowSpec::mToKeep` needs to be updated with the actual index positions in the output table. This was easily fixed by passing in the relevant mapping by reference instead of by const reference. Again, this is a case where I did not see this problem initially because I was dealing solely with MC data. Thanks to Sandro for providing the instructions for how to run a full MC->AOD chain. Also, to reproduce results from `dev`, I had to implement a (faulty) track selection into `AODMcProducerHelpers::updateParticles`. It is important to keep in mind that `AODProducerWorkflowSpec::mToKeep[source][event]` is a mapping from particle number to storage flag. A zero or negative storage flag means that the track is stored. In `dev`, the algorithm is if particle from EG: store it else if particle is physical primary: store it else if particle is physics: store it if particle found in mapping: mark mothers and daughters for storage The important part is the last thing: `if particle found in mapping`. The particle _may_ be stored in the mapping with a negative flag - meaning it should not be stored - which means that we may end up triggering storage of mothers and daughters of a particle that isn't it self stored. In my test of 100 pp events with roughly 100k particles stored in total, this happend 25 times. The correct algorithm is if particle not previously marked for storage and particle is not from EG and particle is not physical primary and particle is not physics: do not store particle go on to next particle store particle and its mothers and daughters In this way, mothers and daughters will _only_ be marked for storage if the particle it self is in fact stored. Currently, the selection implements the `dev` algorithm only so that we can test that `dev` and this MR gives the same results. Once this MR is accepted, the select upgraded to correct algorithm. --- Detectors/AOD/CMakeLists.txt | 4 +- .../AODMcProducerHelpers.h | 324 +++++++++++++ .../AODMcProducerWorkflowSpec.h | 97 +++- .../AODProducerWorkflow/AODProducerHelpers.h | 10 + .../AODProducerWorkflowSpec.h | 91 +++- Detectors/AOD/src/AODMcProducerHelpers.cxx | 431 ++++++++++++++++++ .../AOD/src/AODMcProducerWorkflowSpec.cxx | 414 ++++++++--------- Detectors/AOD/src/AODProducerWorkflowSpec.cxx | 261 +++++------ .../AOD/src/aod-mc-producer-workflow.cxx | 3 + Generators/include/Generators/AODToHepMC.h | 19 +- Generators/src/AODToHepMC.cxx | 40 +- run/CMakeLists.txt | 5 +- run/o2aod_mc_to_hepmc.cxx | 137 +++++- run/o2sim_mctracks_to_aod.cxx | 277 +++++------ 14 files changed, 1530 insertions(+), 583 deletions(-) create mode 100644 Detectors/AOD/include/AODProducerWorkflow/AODMcProducerHelpers.h create mode 100644 Detectors/AOD/src/AODMcProducerHelpers.cxx 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..76a550fe0b893 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,28 @@ 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 +125,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 +150,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 +161,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 +174,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 +185,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 +204,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 +236,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 +// From 636ef7a20b7c30d3dcec8b67ed205d6cb06c4668 Mon Sep 17 00:00:00 2001 From: Christian Holm Christensen Date: Fri, 19 Jan 2024 15:31:50 +0100 Subject: [PATCH 2/2] Code check fix --- Generators/src/AODToHepMC.cxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Generators/src/AODToHepMC.cxx b/Generators/src/AODToHepMC.cxx index 76a550fe0b893..4524bf7f2a60e 100644 --- a/Generators/src/AODToHepMC.cxx +++ b/Generators/src/AODToHepMC.cxx @@ -75,8 +75,9 @@ void AODToHepMC::process(Header const& collision, void AODToHepMC::endEvent() { LOG(debug) << "<<< an event"; - if (not mWriter) + if (not mWriter) { return; + } // If we have a writer, then dump event to output file mWriter->write_event(mEvent); }