diff --git a/Common/Constants/include/CommonConstants/MathConstants.h b/Common/Constants/include/CommonConstants/MathConstants.h index 6870b8ddd5712..9ef3b4dba5ae0 100644 --- a/Common/Constants/include/CommonConstants/MathConstants.h +++ b/Common/Constants/include/CommonConstants/MathConstants.h @@ -22,7 +22,8 @@ namespace constants { namespace math { -constexpr float Almost0 = 1.175494351e-38f; +constexpr float Almost0 = 0x1.0p-126f; // smallest non-denormal float +constexpr float Epsilon = 0x0.000002p0f; // smallest float such that 1 != 1 + Epsilon constexpr float Almost1 = 1.f - 1.0e-6f; constexpr float VeryBig = 1.f / Almost0; diff --git a/Common/Utils/include/CommonUtils/ConfigurableParam.h b/Common/Utils/include/CommonUtils/ConfigurableParam.h index 717a4c425fc82..f44d9efcaea76 100644 --- a/Common/Utils/include/CommonUtils/ConfigurableParam.h +++ b/Common/Utils/include/CommonUtils/ConfigurableParam.h @@ -321,17 +321,19 @@ class ConfigurableParam } // end namespace o2 // a helper macro for boilerplate code in parameter classes -#define O2ParamDef(classname, key) \ - public: \ - classname(TRootIOCtor*) {} \ - classname(classname const&) = delete; \ - \ - private: \ - static constexpr char const* const sKey = key; \ - static classname sInstance; \ - classname() = default; \ - template \ - friend class o2::conf::ConfigurableParamHelper; +#define O2ParamDef(classname, key) \ + public: \ + classname(TRootIOCtor*) {} \ + classname(classname const&) = delete; \ + \ + private: \ + static constexpr char const* const sKey = key; \ + static classname sInstance; \ + classname() = default; \ + template \ + friend class o2::conf::ConfigurableParamHelper; \ + template \ + friend class o2::conf::ConfigurableParamPromoter; // a helper macro to implement necessary symbols in source #define O2ParamImpl(classname) classname classname::sInstance; diff --git a/Common/Utils/include/CommonUtils/ConfigurableParamHelper.h b/Common/Utils/include/CommonUtils/ConfigurableParamHelper.h index 1dc5d5c4c38f8..7d9cb78bb9968 100644 --- a/Common/Utils/include/CommonUtils/ConfigurableParamHelper.h +++ b/Common/Utils/include/CommonUtils/ConfigurableParamHelper.h @@ -45,18 +45,18 @@ class _ParamHelper { private: static std::vector* getDataMembersImpl(std::string const& mainkey, TClass* cl, void*, - std::map const* provmap); + std::map const* provmap, size_t virtualoffset); static void fillKeyValuesImpl(std::string const& mainkey, TClass* cl, void*, boost::property_tree::ptree*, std::map>*, - EnumRegistry*); + EnumRegistry*, size_t offset); static void printWarning(std::type_info const&); static void assignmentImpl(std::string const& mainkey, TClass* cl, void* to, void* from, - std::map* provmap); + std::map* provmap, size_t offset); static void syncCCDBandRegistry(std::string const& mainkey, TClass* cl, void* to, void* from, - std::map* provmap); + std::map* provmap, size_t offset); static void outputMembersImpl(std::ostream& out, std::string const& mainkey, std::vector const* members, bool showProv, bool useLogger); static void printMembersImpl(std::string const& mainkey, std::vector const* members, bool showProv, bool useLogger); @@ -65,6 +65,9 @@ class _ParamHelper template friend class ConfigurableParamHelper; + + template + friend class ConfigurableParamPromoter; }; // ---------------------------------------------------------------- @@ -140,7 +143,7 @@ class ConfigurableParamHelper : virtual public ConfigurableParam return nullptr; } - return _ParamHelper::getDataMembersImpl(getName(), cl, (void*)this, sValueProvenanceMap); + return _ParamHelper::getDataMembersImpl(getName(), cl, (void*)this, sValueProvenanceMap, 0); } // ---------------------------------------------------------------- @@ -153,7 +156,7 @@ class ConfigurableParamHelper : virtual public ConfigurableParam _ParamHelper::printWarning(typeid(P)); return; } - _ParamHelper::fillKeyValuesImpl(getName(), cl, (void*)this, tree, sKeyToStorageMap, sEnumRegistry); + _ParamHelper::fillKeyValuesImpl(getName(), cl, (void*)this, tree, sKeyToStorageMap, sEnumRegistry, 0); } // ---------------------------------------------------------------- @@ -167,7 +170,7 @@ class ConfigurableParamHelper : virtual public ConfigurableParam file->GetObject(getName().c_str(), readback); if (readback != nullptr) { _ParamHelper::assignmentImpl(getName(), TClass::GetClass(typeid(P)), (void*)this, (void*)readback, - sValueProvenanceMap); + sValueProvenanceMap, 0); delete readback; } setRegisterMode(true); @@ -185,7 +188,146 @@ class ConfigurableParamHelper : virtual public ConfigurableParam // setRegisterMode(false); _ParamHelper::syncCCDBandRegistry(getName(), TClass::GetClass(typeid(P)), (void*)this, (void*)externalobj, - sValueProvenanceMap); + sValueProvenanceMap, 0); + setRegisterMode(true); + } + + // ---------------------------------------------------------------- + + void serializeTo(TFile* file) const final + { + file->WriteObjectAny((void*)this, TClass::GetClass(typeid(P)), getName().c_str()); + } +}; + +// Promotes a simple struct Base to a configurable parameter class +// Aka implements all interfaces for a ConfigurableParam P, which shares or +// takes the fields from a Base struct +template +class ConfigurableParamPromoter : public Base, virtual public ConfigurableParam +{ + public: + using ConfigurableParam::ConfigurableParam; + + static const P& Instance() + { + return P::sInstance; + } + + // extracts a copy of the underlying data struct + Base detach() const + { + static_assert(std::copyable, "Base type must be copyable."); + return static_cast(*this); + } + + // ---------------------------------------------------------------- + std::string getName() const final + { + return P::sKey; + } + + // ---------------------------------------------------------------- + // get the provenace of the member with given key + EParamProvenance getMemberProvenance(const std::string& key) const final + { + return getProvenance(getName() + '.' + key); + } + + // ---------------------------------------------------------------- + + // one of the key methods, using introspection to print itself + void printKeyValues(bool showProv = true, bool useLogger = false) const final + { + if (!isInitialized()) { + initialize(); + } + auto members = getDataMembers(); + _ParamHelper::printMembersImpl(getName(), members, showProv, useLogger); + } + + // + size_t getHash() const final + { + return _ParamHelper::getHashImpl(getName(), getDataMembers()); + } + + // ---------------------------------------------------------------- + + void output(std::ostream& out) const final + { + auto members = getDataMembers(); + _ParamHelper::outputMembersImpl(out, getName(), members, true, false); + } + + // ---------------------------------------------------------------- + + // Grab the list of ConfigurableParam data members + // Returns a nullptr if the TClass of the P template class cannot be created. + std::vector* getDataMembers() const + { + // just a helper line to make sure P::sInstance is looked-up + // and that compiler complains about missing static sInstance of type P + // volatile void* ptr = (void*)&P::sInstance; + // static assert on type of sInstance: + static_assert(std::is_same::value, + "static instance must of same type as class"); + + // obtain the TClass for the Base type and delegate further + auto cl = TClass::GetClass(typeid(Base)); + if (!cl) { + _ParamHelper::printWarning(typeid(Base)); + return nullptr; + } + + // we need to put an offset of 8 bytes since internally this is using data members of the Base class + // which doesn't account for the virtual table of P + return _ParamHelper::getDataMembersImpl(getName(), cl, (void*)this, sValueProvenanceMap, 8); + } + + // ---------------------------------------------------------------- + + // fills the data structures with the initial default values + void putKeyValues(boost::property_tree::ptree* tree) final + { + auto cl = TClass::GetClass(typeid(Base)); + if (!cl) { + _ParamHelper::printWarning(typeid(Base)); + return; + } + _ParamHelper::fillKeyValuesImpl(getName(), cl, (void*)this, tree, sKeyToStorageMap, sEnumRegistry, 8); + } + + // ---------------------------------------------------------------- + + void initFrom(TFile* file) final + { + // switch off auto registering since the readback object is + // only a "temporary" singleton + setRegisterMode(false); + P* readback = nullptr; + file->GetObject(getName().c_str(), readback); + if (readback != nullptr) { + _ParamHelper::assignmentImpl(getName(), TClass::GetClass(typeid(Base)), (void*)this, (void*)readback, + sValueProvenanceMap, 8); + delete readback; + } + setRegisterMode(true); + } + + // ---------------------------------------------------------------- + + void syncCCDBandRegistry(void* externalobj) final + { + // We may be getting an external copy from CCDB which is passed as externalobj. + // The task of this function is to + // a) update the internal registry with fields coming from CCDB + // but only if keys have not been modified via RT == command line / ini file + // b) update the external object with with fields having RT provenance + // + setRegisterMode(false); + _ParamHelper::syncCCDBandRegistry(getName(), TClass::GetClass(typeid(Base)), (void*)this, (void*)externalobj, + sValueProvenanceMap, 8); setRegisterMode(true); } diff --git a/Common/Utils/include/CommonUtils/EnumBitOperators.h b/Common/Utils/include/CommonUtils/EnumBitOperators.h new file mode 100644 index 0000000000000..3369a8eacf615 --- /dev/null +++ b/Common/Utils/include/CommonUtils/EnumBitOperators.h @@ -0,0 +1,66 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +#ifndef O2_FRAMEWORK_ENUM_BIT_OPERATORS_H_ +#define O2_FRAMEWORK_ENUM_BIT_OPERATORS_H_ + +#include + +#define O2_DEFINE_ENUM_BIT_OPERATORS(enum_t) \ + constexpr auto operator|(enum_t lhs, enum_t rhs) \ + { \ + return static_cast( \ + static_cast>(lhs) | \ + static_cast>(rhs)); \ + } \ + \ + constexpr auto operator&(enum_t lhs, enum_t rhs) \ + { \ + return static_cast( \ + static_cast>(lhs) & \ + static_cast>(rhs)); \ + } \ + \ + constexpr auto operator^(enum_t lhs, enum_t rhs) \ + { \ + return static_cast( \ + static_cast>(lhs) ^ \ + static_cast>(rhs)); \ + } \ + \ + constexpr auto operator~(enum_t op) \ + { \ + return static_cast( \ + ~static_cast>(op)); \ + } \ + \ + constexpr auto& operator|=(enum_t& lhs, enum_t rhs) \ + { \ + lhs = lhs | rhs; \ + return lhs; \ + } \ + \ + constexpr auto& operator&=(enum_t& lhs, enum_t rhs) \ + { \ + lhs = lhs & rhs; \ + return lhs; \ + } \ + \ + constexpr enum_t& operator^=(enum_t& lhs, enum_t rhs) \ + { \ + lhs = lhs ^ rhs; \ + return lhs; \ + } + +#define O2_ENUM_TEST_BIT(mask, value) ((mask & value) == value) +#define O2_ENUM_SET_BIT(bit) ((1 << bit)) +#define O2_ENUM_ANY_BIT(enum) ((static_cast>(enum) != 0)) + +#endif diff --git a/Common/Utils/include/CommonUtils/TreeStream.h b/Common/Utils/include/CommonUtils/TreeStream.h index 2c55f48c98d3a..2aa02f6509d2c 100644 --- a/Common/Utils/include/CommonUtils/TreeStream.h +++ b/Common/Utils/include/CommonUtils/TreeStream.h @@ -63,6 +63,7 @@ class TreeStream const char* getName() const { return mTree.GetName(); } void setID(int id) { mID = id; } int getID() const { return mID; } + TreeStream& operator<<(const Bool_t& b) { CheckIn('B', &b); @@ -75,6 +76,12 @@ class TreeStream return *this; } + TreeStream& operator<<(const int8_t& i) + { + CheckIn('B', &i); + return *this; + } + TreeStream& operator<<(const UChar_t& c) { CheckIn('b', &c); diff --git a/Common/Utils/src/ConfigurableParamHelper.cxx b/Common/Utils/src/ConfigurableParamHelper.cxx index 0fb213b722e26..f217d402bcb45 100644 --- a/Common/Utils/src/ConfigurableParamHelper.cxx +++ b/Common/Utils/src/ConfigurableParamHelper.cxx @@ -182,19 +182,19 @@ std::string asString(TDataMember const& dm, char* pointer) // potentially other cases to be added here LOG(error) << "COULD NOT REPRESENT AS STRING"; - return nullptr; + return std::string(); } // ---------------------------------------------------------------------- std::vector* _ParamHelper::getDataMembersImpl(std::string const& mainkey, TClass* cl, void* obj, - std::map const* provmap) + std::map const* provmap, size_t globaloffset) { std::vector* members = new std::vector; - auto toDataMember = [&members, obj, mainkey, provmap](const TDataMember* dm, int index, int size) { + auto toDataMember = [&members, obj, mainkey, provmap, globaloffset](const TDataMember* dm, int index, int size) { auto TS = getSizeOfUnderlyingType(*dm); - char* pointer = ((char*)obj) + dm->GetOffset() + index * TS; + char* pointer = ((char*)obj) + dm->GetOffset() + index * TS + globaloffset; const std::string name = getName(dm, index, size); auto value = asString(*dm, pointer); @@ -280,14 +280,14 @@ std::type_info const& nameToTypeInfo(const char* tname, TDataType const* dt) void _ParamHelper::fillKeyValuesImpl(std::string const& mainkey, TClass* cl, void* obj, boost::property_tree::ptree* tree, std::map>* keytostoragemap, - EnumRegistry* enumRegistry) + EnumRegistry* enumRegistry, size_t globaloffset) { boost::property_tree::ptree localtree; - auto fillMap = [obj, &mainkey, &localtree, &keytostoragemap, &enumRegistry](const TDataMember* dm, int index, int size) { + auto fillMap = [obj, &mainkey, &localtree, &keytostoragemap, &enumRegistry, globaloffset](const TDataMember* dm, int index, int size) { const auto name = getName(dm, index, size); auto dt = dm->GetDataType(); auto TS = getSizeOfUnderlyingType(*dm); - char* pointer = ((char*)obj) + dm->GetOffset() + index * TS; + char* pointer = ((char*)obj) + dm->GetOffset() + index * TS + globaloffset; localtree.put(name, asString(*dm, pointer)); auto key = mainkey + "." + name; @@ -355,14 +355,14 @@ bool isMemblockDifferent(char const* block1, char const* block2, int sizeinbytes // ---------------------------------------------------------------------- void _ParamHelper::assignmentImpl(std::string const& mainkey, TClass* cl, void* to, void* from, - std::map* provmap) + std::map* provmap, size_t globaloffset) { - auto assignifchanged = [to, from, &mainkey, provmap](const TDataMember* dm, int index, int size) { + auto assignifchanged = [to, from, &mainkey, provmap, globaloffset](const TDataMember* dm, int index, int size) { const auto name = getName(dm, index, size); auto dt = dm->GetDataType(); auto TS = getSizeOfUnderlyingType(*dm); - char* pointerto = ((char*)to) + dm->GetOffset() + index * TS; - char* pointerfrom = ((char*)from) + dm->GetOffset() + index * TS; + char* pointerto = ((char*)to) + dm->GetOffset() + index * TS + globaloffset; + char* pointerfrom = ((char*)from) + dm->GetOffset() + index * TS + globaloffset; // lambda to update the provenance auto updateProv = [&mainkey, name, provmap]() { @@ -402,14 +402,14 @@ void _ParamHelper::assignmentImpl(std::string const& mainkey, TClass* cl, void* // ---------------------------------------------------------------------- void _ParamHelper::syncCCDBandRegistry(const std::string& mainkey, TClass* cl, void* to, void* from, - std::map* provmap) + std::map* provmap, size_t globaloffset) { - auto sync = [to, from, &mainkey, provmap](const TDataMember* dm, int index, int size) { + auto sync = [to, from, &mainkey, provmap, globaloffset](const TDataMember* dm, int index, int size) { const auto name = getName(dm, index, size); auto dt = dm->GetDataType(); auto TS = getSizeOfUnderlyingType(*dm); - char* pointerto = ((char*)to) + dm->GetOffset() + index * TS; - char* pointerfrom = ((char*)from) + dm->GetOffset() + index * TS; + char* pointerto = ((char*)to) + dm->GetOffset() + index * TS + globaloffset; + char* pointerfrom = ((char*)from) + dm->GetOffset() + index * TS + globaloffset; // check current provenance auto key = mainkey + "." + name; diff --git a/DataFormats/Detectors/CTP/include/DataFormatsCTP/TriggerOffsetsParam.h b/DataFormats/Detectors/CTP/include/DataFormatsCTP/TriggerOffsetsParam.h index f931e9eaa8360..063336e5461ce 100644 --- a/DataFormats/Detectors/CTP/include/DataFormatsCTP/TriggerOffsetsParam.h +++ b/DataFormats/Detectors/CTP/include/DataFormatsCTP/TriggerOffsetsParam.h @@ -24,9 +24,10 @@ namespace ctp struct TriggerOffsetsParam : public o2::conf::ConfigurableParamHelper { static constexpr int MaxNDet = 32; // take with margin to account for possible changes / upgrades int64_t LM_L0 = 15; - int64_t L0_L1 = 280; + int64_t L0_L1 = 281; // trigger input latency int64_t globalInputsShift = 0; // Global shift of inps; customOffset[CTP] is global shift of classes int64_t customOffset[MaxNDet] = {}; + int64_t L0_L1_classes = 280; // trigger input latency O2ParamDef(TriggerOffsetsParam, "TriggerOffsetsParam"); // boilerplate stuff + make principal key }; } // namespace ctp diff --git a/DataFormats/Detectors/TPC/include/DataFormatsTPC/CompressedClusters.h b/DataFormats/Detectors/TPC/include/DataFormatsTPC/CompressedClusters.h index 9f49884035b7e..46da2da2a702e 100644 --- a/DataFormats/Detectors/TPC/include/DataFormatsTPC/CompressedClusters.h +++ b/DataFormats/Detectors/TPC/include/DataFormatsTPC/CompressedClusters.h @@ -29,7 +29,7 @@ struct CompressedClustersCounters { unsigned int nUnattachedClusters = 0; unsigned int nAttachedClustersReduced = 0; unsigned int nSliceRows = 36 * 152; - unsigned char nComppressionModes = 0; + unsigned char nComppressionModes = 0; // Don't fix this name due to ROOT dictionaries! float solenoidBz = -1e6f; int maxTimeBin = -1e6; diff --git a/DataFormats/QualityControl/README.md b/DataFormats/QualityControl/README.md index 486856c983306..33821319b7316 100644 --- a/DataFormats/QualityControl/README.md +++ b/DataFormats/QualityControl/README.md @@ -15,7 +15,7 @@ Data quality is determined through two methods: Both methods utilize the same data format for Flags. During processing (both synchronous and asynchronous), Checks produce Qualities and associate them with Flags. -The Quality Control framework then transmits these Flags to the RCT through a gRPC interface (**not ready yet**, to be done in the scope of QC-978). +The Quality Control framework then transmits these Flags to the RCT through a gRPC interface (more details in QC repository documentation). Detector experts can then review the automatically generated Flags and make any necessary modifications or additions directly in the RCT. ### Quality Control Flag Structure @@ -49,12 +49,13 @@ Each Flag Type has the following attributes: #### Creating and Managing Flag Types * **FlagTypeFactory** ensures a centralized and consistent list of available Flag Types. - New types can only be created through this factory. + New Flags can only be created through this factory. * **[flagTypes.csv](etc/flagTypes.csv)** defines the existing Flag Types, including their ID, name, and "bad quality" determinant, factory method name and a switch to deprecate a flag. The table serves as the source to automatically generate the corresponding methods in FlagTypeFactory. * **Adding new Flag Types:** If a new issue requires a flag not currently defined, propose the addition by contacting the async QC coordinators. They have the authority to add new Flag Types to the RCT. These changes will then be reflected in the [flagTypes.csv](etc/flagTypes.csv) file through a pull request. + Any proposals for new Flag Types should describe the effects on usability of data from analyzer point of view and they should not be detector-specific unless well-argumented. * **Modification of existing Flag Types:** Existing Flag Types should not be modified in terms of their definition. Instead, one may create a new Flag Type and mark the existing one as obsolete in the CSV table. This will add the `[[ deprecated ]]` attribute to the corresponding method. diff --git a/DataFormats/Reconstruction/src/TrackParametrizationWithError.cxx b/DataFormats/Reconstruction/src/TrackParametrizationWithError.cxx index 0dd4a4441c0b3..81963adf79938 100644 --- a/DataFormats/Reconstruction/src/TrackParametrizationWithError.cxx +++ b/DataFormats/Reconstruction/src/TrackParametrizationWithError.cxx @@ -12,14 +12,9 @@ #include "ReconstructionDataFormats/TrackParametrizationWithError.h" #include "ReconstructionDataFormats/Vertex.h" #include "ReconstructionDataFormats/DCA.h" +#include "CommonConstants/MathConstants.h" #include -#ifndef __OPENCL__ -#include -#else -#include -#endif - #ifndef GPUCA_GPUCODE_DEVICE #include #endif @@ -794,11 +789,11 @@ GPUd() auto TrackParametrizationWithError::getPredictedChi2(const Track // get chi2 wrt other track, which must be defined at the same parameters X,alpha // Supplied non-initialized covToSet matrix is filled by inverse combined matrix for further use - if (gpu::CAMath::Abs(this->getAlpha() - rhs.getAlpha()) > FLT_EPSILON) { + if (gpu::CAMath::Abs(this->getAlpha() - rhs.getAlpha()) > o2::constants::math::Epsilon) { LOG(error) << "The reference Alpha of the tracks differ: " << this->getAlpha() << " : " << rhs.getAlpha(); return 2.f * HugeF; } - if (gpu::CAMath::Abs(this->getX() - rhs.getX()) > FLT_EPSILON) { + if (gpu::CAMath::Abs(this->getX() - rhs.getX()) > o2::constants::math::Epsilon) { LOG(error) << "The reference X of the tracks differ: " << this->getX() << " : " << rhs.getX(); return 2.f * HugeF; } @@ -827,11 +822,11 @@ GPUd() bool TrackParametrizationWithError::update(const TrackParametriz // update track with other track, the inverted combined cov matrix should be supplied // consider skipping this check, since it is usually already done upstream - if (gpu::CAMath::Abs(this->getAlpha() - rhs.getAlpha()) > FLT_EPSILON) { + if (gpu::CAMath::Abs(this->getAlpha() - rhs.getAlpha()) > o2::constants::math::Epsilon) { LOG(error) << "The reference Alpha of the tracks differ: " << this->getAlpha() << " : " << rhs.getAlpha(); return false; } - if (gpu::CAMath::Abs(this->getX() - rhs.getX()) > FLT_EPSILON) { + if (gpu::CAMath::Abs(this->getX() - rhs.getX()) > o2::constants::math::Epsilon) { LOG(error) << "The reference X of the tracks differ: " << this->getX() << " : " << rhs.getX(); return false; } diff --git a/DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h b/DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h index 229a1a7a8a535..6b1690946e951 100644 --- a/DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h +++ b/DataFormats/simulation/include/SimulationDataFormat/O2DatabasePDG.h @@ -235,62 +235,38 @@ inline void O2DatabasePDG::addALICEParticles(TDatabasePDG* db) //Hyper nuclei and exotica ionCode = 1010010030; if (!db->GetParticle(ionCode)) { - db->AddParticle("HyperTriton", "HyperTriton", 2.99131, kFALSE, + db->AddParticle("HyperTriton", "HyperTriton", 2.991134, kFALSE, 2.5e-15, 3, "Ion", ionCode); } ionCode = -1010010030; if (!db->GetParticle(ionCode)) { - db->AddParticle("AntiHyperTriton", "AntiHyperTriton", 2.99131, kFALSE, + db->AddParticle("AntiHyperTriton", "AntiHyperTriton", 2.991134, kFALSE, 2.5e-15, 3, "Ion", ionCode); } //hyper hydrogen 4 ground state ionCode = 1010010040; if (!db->GetParticle(ionCode)) { - db->AddParticle("Hyperhydrog4", "Hyperhydrog4", 3.9226, kFALSE, + db->AddParticle("Hyperhydrog4", "Hyperhydrog4", 3.922434, kFALSE, 2.5e-15, 3, "Ion", ionCode); } //anti hyper hydrogen 4 ground state ionCode = -1010010040; if (!db->GetParticle(ionCode)) { - db->AddParticle("AntiHyperhydrog4", "AntiHyperhydrog4", 3.9226, kFALSE, - 2.5e-15, 3, "Ion", ionCode); - } - //hyper hydrogen 4 excited state - ionCode = 1010010041; - if (!db->GetParticle(ionCode)) { - db->AddParticle("Hyperhydrog4*", "Hyperhydrog4*", 3.9237, kFALSE, - 2.5e-15, 3, "Ion", ionCode); - } - //anti hyper hydrogen 4 excited state - ionCode = -1010010041; - if (!db->GetParticle(ionCode)) { - db->AddParticle("AntiHyperhydrog4*", "AntiHyperhydrog4*", 3.9237, kFALSE, + db->AddParticle("AntiHyperhydrog4", "AntiHyperhydrog4", 3.922434, kFALSE, 2.5e-15, 3, "Ion", ionCode); } //hyper helium 4 ground state ionCode = 1010020040; if (!db->GetParticle(ionCode)) { - db->AddParticle("Hyperhelium4", "Hyperhelium4", 3.9217, kFALSE, + db->AddParticle("Hyperhelium4", "Hyperhelium4", 3.921728, kFALSE, 2.5e-15, 6, "Ion", ionCode); } //anti hyper helium 4 ground state ionCode = -1010020040; if (!db->GetParticle(ionCode)) { - db->AddParticle("AntiHyperhelium4", "AntiHyperhelium4", 3.9217, kFALSE, - 2.5e-15, 6, "Ion", ionCode); - } - //hyper helium 4 excited state - ionCode = 1010020041; - if (!db->GetParticle(ionCode)) { - db->AddParticle("Hyperhelium4*", "Hyperhelium4*", 3.9231, kFALSE, - 2.5e-15, 6, "Ion", ionCode); - } - //anti hyper helium 4 excited state - ionCode = -1010020041; - if (!db->GetParticle(ionCode)) { - db->AddParticle("AntiHyperhelium4*", "AntiHyperhelium4*", 3.9231, kFALSE, + db->AddParticle("AntiHyperhelium4", "AntiHyperhelium4", 3.921728, kFALSE, 2.5e-15, 6, "Ion", ionCode); } @@ -309,13 +285,13 @@ inline void O2DatabasePDG::addALICEParticles(TDatabasePDG* db) ionCode = 1010020050; if (!db->GetParticle(ionCode)) { - db->AddParticle("Hyperhelium5", "Hyperhelium5", 4.841, kFALSE, + db->AddParticle("Hyperhelium5", "Hyperhelium5", 4.839961, kFALSE, 2.5e-15, 6, "Ion", ionCode); } ionCode = -1010020050; if (!db->GetParticle(ionCode)) { - db->AddParticle("AntiHyperhelium5", "AntiHyperhelium5", 4.841, kFALSE, + db->AddParticle("AntiHyperhelium5", "AntiHyperhelium5", 4.839961, kFALSE, 2.5e-15, 6, "Ion", ionCode); } @@ -331,6 +307,20 @@ inline void O2DatabasePDG::addALICEParticles(TDatabasePDG* db) 2.5e-15, 6, "Ion", ionCode); } + // 4-Xi-He + ionCode = 1120020040; + if (!db->GetParticle(ionCode)) { + db->AddParticle("4XiHe", "4XiHe", 4.128, kFALSE, 4.04e-15, 3, "Ion", ionCode); + db->AddAntiParticle("Anti4XiHe", -ionCode); + } + + // 4-Xi-H + ionCode = 1120010040; + if (!db->GetParticle(ionCode)) { + db->AddParticle("4XiH", "4XiH", 4.128, kFALSE, 4.04e-15, 3, "Ion", ionCode); + db->AddAntiParticle("Anti4XiH", -ionCode); + } + // hyper helium 4 sigma ionCode = 1110020040; if (!db->GetParticle(ionCode)) { diff --git a/Detectors/AOD/include/AODProducerWorkflow/AODProducerWorkflowSpec.h b/Detectors/AOD/include/AODProducerWorkflow/AODProducerWorkflowSpec.h index 94f4526fe30a1..d9481917f9a05 100644 --- a/Detectors/AOD/include/AODProducerWorkflow/AODProducerWorkflowSpec.h +++ b/Detectors/AOD/include/AODProducerWorkflow/AODProducerWorkflowSpec.h @@ -30,7 +30,11 @@ #include "TStopwatch.h" #include "ZDCBase/Constants.h" #include "GlobalTracking/MatchGlobalFwd.h" +#include "CommonUtils/TreeStreamRedirector.h" +#include "CommonUtils/EnumBitOperators.h" +#include +#include #include #include #include @@ -203,7 +207,15 @@ class BunchCrossings std::vector mTimeWindows; // the time window structure covering the complete duration of mBCTimeVector double mWindowSize; // the size of a single time window -}; // end internal class +}; // end internal class + +// Steering bits for additional output during AOD production +enum struct AODProducerStreamerMask : uint8_t { + None = 0, + TrackQA = O2_ENUM_SET_BIT(0), + All = std::numeric_limits>::max(), +}; +O2_DEFINE_ENUM_BIT_OPERATORS(AODProducerStreamerMask) class AODProducerWorkflowDPL : public Task { @@ -241,6 +253,9 @@ class AODProducerWorkflowDPL : public Task std::unordered_set mGIDUsedBySVtx; std::unordered_set mGIDUsedByStr; + AODProducerStreamerMask mStreamerMask; + std::shared_ptr mStreamer; + int mNThreads = 1; bool mUseMC = true; bool mEnableSV = true; // enable secondary vertices @@ -339,6 +354,7 @@ class AODProducerWorkflowDPL : public Task uint32_t mTrackCovOffDiag = 0xFFFF0000; // 7 bits uint32_t mTrackSignal = 0xFFFFFF00; // 15 bits uint32_t mTrackTime = 0xFFFFFFFF; // use full float precision for time + uint32_t mTPCTime0 = 0xFFFFFFE0; // 18 bits, providing 14256./(1<<19) = 0.027 TB precision e.g., ~0.13 mm in z uint32_t mTrackTimeError = 0xFFFFFF00; // 15 bits uint32_t mTrackPosEMCAL = 0xFFFFFF00; // 15 bits uint32_t mTracklets = 0xFFFFFF00; // 15 bits @@ -397,18 +413,28 @@ class AODProducerWorkflowDPL : public Task struct TrackQA { GID trackID; - float tpcTime0; - int16_t tpcdcaR; - int16_t tpcdcaZ; - uint8_t tpcClusterByteMask; - uint8_t tpcdEdxMax0R; - uint8_t tpcdEdxMax1R; - uint8_t tpcdEdxMax2R; - uint8_t tpcdEdxMax3R; - uint8_t tpcdEdxTot0R; - uint8_t tpcdEdxTot1R; - uint8_t tpcdEdxTot2R; - uint8_t tpcdEdxTot3R; + float tpcTime0{}; + int16_t tpcdcaR{}; + int16_t tpcdcaZ{}; + uint8_t tpcClusterByteMask{}; + uint8_t tpcdEdxMax0R{}; + uint8_t tpcdEdxMax1R{}; + uint8_t tpcdEdxMax2R{}; + uint8_t tpcdEdxMax3R{}; + uint8_t tpcdEdxTot0R{}; + uint8_t tpcdEdxTot1R{}; + uint8_t tpcdEdxTot2R{}; + uint8_t tpcdEdxTot3R{}; + int8_t dRefContY{std::numeric_limits::min()}; + int8_t dRefContZ{std::numeric_limits::min()}; + int8_t dRefContSnp{std::numeric_limits::min()}; + int8_t dRefContTgl{std::numeric_limits::min()}; + int8_t dRefContQ2Pt{std::numeric_limits::min()}; + int8_t dRefGloY{std::numeric_limits::min()}; + int8_t dRefGloZ{std::numeric_limits::min()}; + int8_t dRefGloSnp{std::numeric_limits::min()}; + int8_t dRefGloTgl{std::numeric_limits::min()}; + int8_t dRefGloQ2Pt{std::numeric_limits::min()}; }; // helper struct for addToFwdTracksTable() diff --git a/Detectors/AOD/src/AODProducerWorkflowSpec.cxx b/Detectors/AOD/src/AODProducerWorkflowSpec.cxx index 6c3a418612478..8a2443b57c7ff 100644 --- a/Detectors/AOD/src/AODProducerWorkflowSpec.cxx +++ b/Detectors/AOD/src/AODProducerWorkflowSpec.cxx @@ -51,6 +51,7 @@ #include "Framework/DataTypes.h" #include "Framework/TableBuilder.h" #include "Framework/CCDBParamSpec.h" +#include "CommonUtils/TreeStreamRedirector.h" #include "FT0Base/Geometry.h" #include "GlobalTracking/MatchTOF.h" #include "ReconstructionDataFormats/Cascade.h" @@ -85,8 +86,10 @@ #include "MathUtils/Utils.h" #include "Math/SMatrix.h" #include "TString.h" +#include #include #include +#include #include #include #include @@ -355,25 +358,47 @@ void AODProducerWorkflowDPL::addToTracksExtraTable(TracksExtraCursorType& tracks template void AODProducerWorkflowDPL::addToTracksQATable(TracksQACursorType& tracksQACursor, TrackQA& trackQAInfoHolder) { - - // trackQA - tracksQACursor( - - // truncateFloatFraction(trackQAInfoHolder.tpcdcaR, mTrackChi2), - // truncateFloatFraction(trackQAInfoHolder.tpcdcaZ, mTrackChi2), - trackQAInfoHolder.trackID, - trackQAInfoHolder.tpcTime0, - trackQAInfoHolder.tpcdcaR, - trackQAInfoHolder.tpcdcaZ, - trackQAInfoHolder.tpcClusterByteMask, - trackQAInfoHolder.tpcdEdxMax0R, - trackQAInfoHolder.tpcdEdxMax1R, - trackQAInfoHolder.tpcdEdxMax2R, - trackQAInfoHolder.tpcdEdxMax3R, - trackQAInfoHolder.tpcdEdxTot0R, - trackQAInfoHolder.tpcdEdxTot1R, - trackQAInfoHolder.tpcdEdxTot2R, - trackQAInfoHolder.tpcdEdxTot3R); + if constexpr (std::is_same_v) { // TODO remove remove once version changes + tracksQACursor( + trackQAInfoHolder.trackID, + truncateFloatFraction(trackQAInfoHolder.tpcTime0, mTPCTime0), + trackQAInfoHolder.tpcdcaR, + trackQAInfoHolder.tpcdcaZ, + trackQAInfoHolder.tpcClusterByteMask, + trackQAInfoHolder.tpcdEdxMax0R, + trackQAInfoHolder.tpcdEdxMax1R, + trackQAInfoHolder.tpcdEdxMax2R, + trackQAInfoHolder.tpcdEdxMax3R, + trackQAInfoHolder.tpcdEdxTot0R, + trackQAInfoHolder.tpcdEdxTot1R, + trackQAInfoHolder.tpcdEdxTot2R, + trackQAInfoHolder.tpcdEdxTot3R, + trackQAInfoHolder.dRefContY, + trackQAInfoHolder.dRefContZ, + trackQAInfoHolder.dRefContSnp, + trackQAInfoHolder.dRefContTgl, + trackQAInfoHolder.dRefContQ2Pt, + trackQAInfoHolder.dRefGloY, + trackQAInfoHolder.dRefGloZ, + trackQAInfoHolder.dRefGloSnp, + trackQAInfoHolder.dRefGloTgl, + trackQAInfoHolder.dRefGloQ2Pt); + } else { + tracksQACursor( + trackQAInfoHolder.trackID, + trackQAInfoHolder.tpcTime0, + trackQAInfoHolder.tpcdcaR, + trackQAInfoHolder.tpcdcaZ, + trackQAInfoHolder.tpcClusterByteMask, + trackQAInfoHolder.tpcdEdxMax0R, + trackQAInfoHolder.tpcdEdxMax1R, + trackQAInfoHolder.tpcdEdxMax2R, + trackQAInfoHolder.tpcdEdxMax3R, + trackQAInfoHolder.tpcdEdxTot0R, + trackQAInfoHolder.tpcdEdxTot1R, + trackQAInfoHolder.tpcdEdxTot2R, + trackQAInfoHolder.tpcdEdxTot3R); + } } template @@ -1664,6 +1689,14 @@ void AODProducerWorkflowDPL::init(InitContext& ic) mThinTracks = ic.options().get("thin-tracks"); mPropTracks = ic.options().get("propagate-tracks"); mPropMuons = ic.options().get("propagate-muons"); + if (auto s = ic.options().get("with-streamers"); !s.empty()) { + mStreamerMask = static_cast(std::stoul(s, nullptr, 2)); + if (O2_ENUM_ANY_BIT(mStreamerMask)) { + LOGP(info, "Writing streamer data with mask {:0{}b}", static_cast>(mStreamerMask), std::numeric_limits>::digits); + } else { + LOGP(warn, "Specified non-default empty streamer mask!"); + } + } mTrackQCFraction = ic.options().get("trackqc-fraction"); mTrackQCNTrCut = ic.options().get("trackqc-NTrCut"); if (auto seed = ic.options().get("seed"); seed == 0) { @@ -1705,6 +1738,7 @@ void AODProducerWorkflowDPL::init(InitContext& ic) mTrackCovOffDiag = 0xFFFFFFFF; mTrackSignal = 0xFFFFFFFF; mTrackTime = 0xFFFFFFFF; + mTPCTime0 = 0xFFFFFFFF; mTrackTimeError = 0xFFFFFFFF; mTrackPosEMCAL = 0xFFFFFFFF; mTracklets = 0xFFFFFFFF; @@ -1748,6 +1782,10 @@ void AODProducerWorkflowDPL::init(InitContext& ic) mHeavyIonUpdate = when; mTimer.Reset(); + + if (O2_ENUM_ANY_BIT(mStreamerMask)) { + mStreamer = std::make_unique("AO2DStreamer.root", "RECREATE"); + } } void AODProducerWorkflowDPL::run(ProcessingContext& pc) @@ -1816,7 +1854,7 @@ void AODProducerWorkflowDPL::run(ProcessingContext& pc) auto tracksCursor = createTableCursor(pc); auto tracksCovCursor = createTableCursor(pc); auto tracksExtraCursor = createTableCursor(pc); - auto tracksQACursor = createTableCursor(pc); + auto tracksQACursor = createTableCursor(pc); auto ambigTracksCursor = createTableCursor(pc); auto ambigMFTTracksCursor = createTableCursor(pc); auto ambigFwdTracksCursor = createTableCursor(pc); @@ -2534,16 +2572,15 @@ AODProducerWorkflowDPL::TrackQA AODProducerWorkflowDPL::processBarrelTrackQA(int TrackQA trackQAHolder; auto contributorsGID = data.getTPCContributorGID(trackIndex); const auto& trackPar = data.getTrackParam(trackIndex); - // auto src = trackIndex.getSource(); if (contributorsGID.isIndexSet()) { + auto prop = o2::base::Propagator::Instance(); const auto& tpcOrig = data.getTPCTrack(contributorsGID); /// getDCA - should be done with the copy of TPC only track - // LOGP(info, "GloIdx: {} TPCIdx: {}, NTPCTracks: {}", trackIndex.asString(), contributorsGID.asString(), data.getTPCTracks().size()); - o2::track::TrackParametrization tpcTMP = tpcOrig; /// get backup of the track - o2::base::Propagator::MatCorrType mMatType = o2::base::Propagator::MatCorrType::USEMatCorrLUT; /// should be parameterized - o2::dataformats::VertexBase v = mVtx.getMeanVertex(collisionID < 0 ? 0.f : data.getPrimaryVertex(collisionID).getZ()); + o2::track::TrackParametrization tpcTMP = tpcOrig; /// get backup of the track + const o2::base::Propagator::MatCorrType mMatType = o2::base::Propagator::MatCorrType::USEMatCorrLUT; /// should be parameterized + const o2::dataformats::VertexBase v = mVtx.getMeanVertex(collisionID < 0 ? 0.f : data.getPrimaryVertex(collisionID).getZ()); o2::gpu::gpustd::array dcaInfo{-999., -999.}; - if (o2::base::Propagator::Instance()->propagateToDCABxByBz({v.getX(), v.getY(), v.getZ()}, tpcTMP, 2.f, mMatType, &dcaInfo)) { + if (prop->propagateToDCABxByBz({v.getX(), v.getY(), v.getZ()}, tpcTMP, 2.f, mMatType, &dcaInfo)) { trackQAHolder.tpcdcaR = 100. * dcaInfo[0] / sqrt(1. + trackPar.getQ2Pt() * trackPar.getQ2Pt()); trackQAHolder.tpcdcaZ = 100. * dcaInfo[1] / sqrt(1. + trackPar.getQ2Pt() * trackPar.getQ2Pt()); } @@ -2567,7 +2604,7 @@ AODProducerWorkflowDPL::TrackQA AODProducerWorkflowDPL::processBarrelTrackQA(int } trackQAHolder.tpcTime0 = tpcOrig.getTime0(); trackQAHolder.tpcClusterByteMask = byteMask; - float dEdxNorm = (tpcOrig.getdEdx().dEdxTotTPC > 0) ? 100. / tpcOrig.getdEdx().dEdxTotTPC : 0; + const float dEdxNorm = (tpcOrig.getdEdx().dEdxTotTPC > 0) ? 100. / tpcOrig.getdEdx().dEdxTotTPC : 0; trackQAHolder.tpcdEdxMax0R = uint8_t(tpcOrig.getdEdx().dEdxMaxIROC * dEdxNorm); trackQAHolder.tpcdEdxMax1R = uint8_t(tpcOrig.getdEdx().dEdxMaxOROC1 * dEdxNorm); trackQAHolder.tpcdEdxMax2R = uint8_t(tpcOrig.getdEdx().dEdxMaxOROC2 * dEdxNorm); @@ -2577,7 +2614,99 @@ AODProducerWorkflowDPL::TrackQA AODProducerWorkflowDPL::processBarrelTrackQA(int trackQAHolder.tpcdEdxTot1R = uint8_t(tpcOrig.getdEdx().dEdxTotOROC1 * dEdxNorm); trackQAHolder.tpcdEdxTot2R = uint8_t(tpcOrig.getdEdx().dEdxTotOROC2 * dEdxNorm); trackQAHolder.tpcdEdxTot3R = uint8_t(tpcOrig.getdEdx().dEdxTotOROC3 * dEdxNorm); - /// + + if constexpr (std::is_same_v) { // TODO remove remove once version changes + // Add matching information at a reference point (defined by + // o2::aod::track::trackQARefRadius) in the same frame as the global track + // without material corrections and error propagation + if (auto itsContGID = data.getITSContributorGID(trackIndex); itsContGID.isIndexSet() && itsContGID.getSource() != GIndex::ITSAB) { + const auto& itsOrig = data.getITSTrack(itsContGID); + o2::track::TrackPar gloCopy = trackPar; + o2::track::TrackPar itsCopy = itsOrig; + o2::track::TrackPar tpcCopy = tpcOrig; + if (prop->propagateToX(gloCopy, o2::aod::track::trackQARefRadius, prop->getNominalBz(), o2::base::Propagator::MAX_SIN_PHI, o2::base::Propagator::MAX_STEP, mMatCorr) && + prop->propagateToAlphaX(tpcCopy, gloCopy.getAlpha(), o2::aod::track::trackQARefRadius, false, o2::base::Propagator::MAX_SIN_PHI, o2::base::Propagator::MAX_STEP, 1, mMatCorr) && + prop->propagateToAlphaX(itsCopy, gloCopy.getAlpha(), o2::aod::track::trackQARefRadius, false, o2::base::Propagator::MAX_SIN_PHI, o2::base::Propagator::MAX_STEP, 1, mMatCorr)) { + // All tracks are now at the same radius and in the same frame and we can calculate the deltas wrt. to the global track + // The scale is defined by the global track scaling depending on beta0 + const float beta0 = std::sqrt(std::min(50.f / tpcOrig.getdEdx().dEdxMaxTPC, 1.f)); + const float qpt = gloCopy.getQ2Pt(); + const float x = qpt / beta0; + // scaling is defined as sigmaBins/sqrt(p0^2 + (p1 * q/pt / beta)^2) + auto scaleCont = [&x](int i) -> float { + return o2::aod::track::trackQAScaleBins / std::sqrt(o2::aod::track::trackQAScaleContP0[i] * o2::aod::track::trackQAScaleContP0[i] + (o2::aod::track::trackQAScaleContP1[i] * x) * (o2::aod::track::trackQAScaleContP1[i] * x)); + }; + auto scaleGlo = [&x](int i) -> float { + return o2::aod::track::trackQAScaleBins / std::sqrt(o2::aod::track::trackQAScaleGloP0[i] * o2::aod::track::trackQAScaleGloP0[i] + (o2::aod::track::trackQAScaleGloP1[i] * x) * (o2::aod::track::trackQAScaleGloP1[i] * x)); + }; + + // This allows to safely clamp any float to one byte, using the + // minmal/maximum values as under-/overflow borders and rounding to the nearest integer + auto safeInt8Clamp = [](auto value) -> int8_t { + using ValType = decltype(value); + return static_cast(TMath::Nint(std::clamp(value, static_cast(std::numeric_limits::min()), static_cast(std::numeric_limits::max())))); + }; + + // Calculate deltas for contributors + trackQAHolder.dRefContY = safeInt8Clamp((itsCopy.getY() - tpcCopy.getY()) * scaleCont(0)); + trackQAHolder.dRefContZ = safeInt8Clamp((itsCopy.getZ() - tpcCopy.getZ()) * scaleCont(1)); + trackQAHolder.dRefContSnp = safeInt8Clamp((itsCopy.getSnp() - tpcCopy.getSnp()) * scaleCont(2)); + trackQAHolder.dRefContTgl = safeInt8Clamp((itsCopy.getTgl() - tpcCopy.getTgl()) * scaleCont(3)); + trackQAHolder.dRefContQ2Pt = safeInt8Clamp((itsCopy.getQ2Pt() - tpcCopy.getQ2Pt()) * scaleCont(4)); + // Calculate deltas for global track against averaged contributors + trackQAHolder.dRefGloY = safeInt8Clamp(((itsCopy.getY() + tpcCopy.getY()) * 0.5f - gloCopy.getY()) * scaleGlo(0)); + trackQAHolder.dRefGloZ = safeInt8Clamp(((itsCopy.getZ() + tpcCopy.getZ()) * 0.5f - gloCopy.getZ()) * scaleGlo(1)); + trackQAHolder.dRefGloSnp = safeInt8Clamp(((itsCopy.getSnp() + tpcCopy.getSnp()) * 0.5f - gloCopy.getSnp()) * scaleGlo(2)); + trackQAHolder.dRefGloTgl = safeInt8Clamp(((itsCopy.getTgl() + tpcCopy.getTgl()) * 0.5f - gloCopy.getTgl()) * scaleGlo(3)); + trackQAHolder.dRefGloQ2Pt = safeInt8Clamp(((itsCopy.getQ2Pt() + tpcCopy.getQ2Pt()) * 0.5f - gloCopy.getQ2Pt()) * scaleGlo(4)); + + if (O2_ENUM_TEST_BIT(mStreamerMask, AODProducerStreamerMask::TrackQA)) { + (*mStreamer) << "trackQA" + << "trackITSOrig=" << itsOrig + << "trackTPCOrig=" << tpcOrig + << "trackITSTPCOrig=" << trackPar + << "trackITSProp=" << itsCopy + << "trackTPCProp=" << tpcCopy + << "trackITSTPCProp=" << gloCopy + << "refRadius=" << o2::aod::track::trackQARefRadius + << "scaleBins=" << o2::aod::track::trackQAScaleBins + << "scaleCont0=" << scaleCont(0) + << "scaleCont1=" << scaleCont(1) + << "scaleCont2=" << scaleCont(2) + << "scaleCont3=" << scaleCont(3) + << "scaleCont4=" << scaleCont(4) + << "scaleGlo0=" << scaleGlo(0) + << "scaleGlo1=" << scaleGlo(1) + << "scaleGlo2=" << scaleGlo(2) + << "scaleGlo3=" << scaleGlo(3) + << "scaleGlo4=" << scaleGlo(4) + << "trackQAHolder.tpcTime0=" << trackQAHolder.tpcTime0 + << "trackQAHolder.tpcdcaR=" << trackQAHolder.tpcdcaR + << "trackQAHolder.tpcdcaZ=" << trackQAHolder.tpcdcaZ + << "trackQAHolder.tpcdcaClusterByteMask=" << trackQAHolder.tpcClusterByteMask + << "trackQAHolder.tpcdEdxMax0R=" << trackQAHolder.tpcdEdxMax0R + << "trackQAHolder.tpcdEdxMax1R=" << trackQAHolder.tpcdEdxMax1R + << "trackQAHolder.tpcdEdxMax2R=" << trackQAHolder.tpcdEdxMax2R + << "trackQAHolder.tpcdEdxMax3R=" << trackQAHolder.tpcdEdxMax3R + << "trackQAHolder.tpcdEdxTot0R=" << trackQAHolder.tpcdEdxTot0R + << "trackQAHolder.tpcdEdxTot1R=" << trackQAHolder.tpcdEdxTot1R + << "trackQAHolder.tpcdEdxTot2R=" << trackQAHolder.tpcdEdxTot2R + << "trackQAHolder.tpcdEdxTot3R=" << trackQAHolder.tpcdEdxTot3R + << "trackQAHolder.dRefContY=" << trackQAHolder.dRefContY + << "trackQAHolder.dRefContZ=" << trackQAHolder.dRefContZ + << "trackQAHolder.dRefContSnp=" << trackQAHolder.dRefContSnp + << "trackQAHolder.dRefContTgl=" << trackQAHolder.dRefContTgl + << "trackQAHolder.dRefContQ2Pt=" << trackQAHolder.dRefContQ2Pt + << "trackQAHolder.dRefGloY=" << trackQAHolder.dRefGloY + << "trackQAHolder.dRefGloZ=" << trackQAHolder.dRefGloZ + << "trackQAHolder.dRefGloSnp=" << trackQAHolder.dRefGloSnp + << "trackQAHolder.dRefGloTgl=" << trackQAHolder.dRefGloTgl + << "trackQAHolder.dRefGloQ2Pt=" << trackQAHolder.dRefGloQ2Pt + << "\n"; + } + } + } + } } return trackQAHolder; @@ -2944,6 +3073,8 @@ void AODProducerWorkflowDPL::endOfStream(EndOfStreamContext& /*ec*/) { LOGF(info, "aod producer dpl total timing: Cpu: %.3e Real: %.3e s in %d slots", mTimer.CpuTime(), mTimer.RealTime(), mTimer.Counter() - 1); + + mStreamer.reset(); } DataProcessorSpec getAODProducerWorkflowSpec(GID::mask_t src, bool enableSV, bool enableStrangenessTracking, bool useMC, bool CTPConfigPerRun) @@ -3076,6 +3207,7 @@ DataProcessorSpec getAODProducerWorkflowSpec(GID::mask_t src, bool enableSV, boo ConfigParamSpec{"thin-tracks", VariantType::Bool, false, {"Produce thinned track tables"}}, 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"}}, + ConfigParamSpec{"with-streamers", VariantType::String, "", {"Bit-mask to steer writing of intermediate streamer files"}}, ConfigParamSpec{"seed", VariantType::Int, 0, {"Set seed for random generator used for sampling (0 (default) means using a random_device)"}}, }}; } diff --git a/Detectors/CTP/reconstruction/include/CTPReconstruction/CTFCoder.h b/Detectors/CTP/reconstruction/include/CTPReconstruction/CTFCoder.h index 6ffb3575207e5..9189df5d12685 100644 --- a/Detectors/CTP/reconstruction/include/CTPReconstruction/CTFCoder.h +++ b/Detectors/CTP/reconstruction/include/CTPReconstruction/CTFCoder.h @@ -25,6 +25,7 @@ #include "DetectorsBase/CTFCoderBase.h" #include "CTPReconstruction/CTFHelper.h" #include "CTPReconstruction/RawDataDecoder.h" +#include "DataFormatsCTP/Configuration.h" class TTree; @@ -53,6 +54,9 @@ class CTFCoder : public o2::ctf::CTFCoderBase void createCoders(const std::vector& bufVec, o2::ctf::CTFCoderBase::OpType op) final; void setDecodeInps(bool decodeinps) { mDecodeInps = decodeinps; } + void setCTPConfig(CTPConfiguration cfg) { mCTPConfig = std::move(cfg); } + bool getDecodeInps() { return mDecodeInps; } + CTPConfiguration& getCTPConfig() { return mCTPConfig; } bool canApplyBCShiftInputs(const o2::InteractionRecord& ir) const { return canApplyBCShift(ir, mBCShiftInputs); } private: @@ -62,6 +66,7 @@ class CTFCoder : public o2::ctf::CTFCoderBase void appendToTree(TTree& tree, CTF& ec); void readFromTree(TTree& tree, int entry, std::vector& data, LumiInfo& lumi); std::vector mDataFilt; + CTPConfiguration mCTPConfig; int mBCShiftInputs = 0; bool mDecodeInps = false; }; @@ -215,8 +220,13 @@ o2::ctf::CTFIOSize CTFCoder::decode(const CTF::base& ec, VTRG& data, LumiInfo& l } } if (mDecodeInps) { + uint64_t trgclassmask = 0xffffffffffffffff; + if (mCTPConfig.getRunNumber() != 0) { + trgclassmask = mCTPConfig.getTriggerClassMask(); + } + // std::cout << "trgclassmask:" << std::hex << trgclassmask << std::dec << std::endl; o2::pmr::vector digits; - o2::ctp::RawDataDecoder::shiftInputs(digitsMap, digits, mFirstTFOrbit); + o2::ctp::RawDataDecoder::shiftInputs(digitsMap, digits, mFirstTFOrbit, trgclassmask); for (auto const& dig : digits) { data.emplace_back(dig); } diff --git a/Detectors/CTP/reconstruction/include/CTPReconstruction/RawDataDecoder.h b/Detectors/CTP/reconstruction/include/CTPReconstruction/RawDataDecoder.h index c50079f9f8717..16a8ec6a6bef1 100644 --- a/Detectors/CTP/reconstruction/include/CTPReconstruction/RawDataDecoder.h +++ b/Detectors/CTP/reconstruction/include/CTPReconstruction/RawDataDecoder.h @@ -22,6 +22,7 @@ #include "Framework/InputRecord.h" #include "DataFormatsCTP/Digits.h" #include "DataFormatsCTP/LumiInfo.h" +#include "DataFormatsCTP/Configuration.h" namespace o2 { @@ -43,14 +44,16 @@ class RawDataDecoder void setVerbose(bool v) { mVerbose = v; } void setMAXErrors(int m) { mErrorMax = m; } int setLumiInp(int lumiinp, std::string inp); + void setCTPConfig(CTPConfiguration cfg) { mCTPConfig = std::move(cfg); }; uint32_t getIRRejected() const { return mIRRejected; } uint32_t getTCRRejected() const { return mTCRRejected; } std::vector& getTFOrbits() { return mTFOrbits; } int getErrorIR() { return mErrorIR; } int getErrorTCR() { return mErrorTCR; } + CTPConfiguration& getCTPConfig() { return mCTPConfig; } int init(); static int shiftNew(const o2::InteractionRecord& irin, uint32_t TFOrbit, std::bitset<48>& inpmask, int64_t shift, int level, std::map& digmap); - static int shiftInputs(std::map& digitsMap, o2::pmr::vector& digits, uint32_t TFOrbit); + static int shiftInputs(std::map& digitsMap, o2::pmr::vector& digits, uint32_t TFOrbit, uint64_t trgclassmask = 0xffffffffffffffff); private: static constexpr uint32_t TF_TRIGGERTYPE_MASK = 0x800; @@ -79,6 +82,7 @@ class RawDataDecoder int mErrorTCR = 0; int mErrorMax = 3; bool mStickyError = false; + CTPConfiguration mCTPConfig; }; } // namespace ctp } // namespace o2 diff --git a/Detectors/CTP/reconstruction/src/RawDataDecoder.cxx b/Detectors/CTP/reconstruction/src/RawDataDecoder.cxx index 4e3d480e463cd..74e5b7481163d 100644 --- a/Detectors/CTP/reconstruction/src/RawDataDecoder.cxx +++ b/Detectors/CTP/reconstruction/src/RawDataDecoder.cxx @@ -89,7 +89,7 @@ int RawDataDecoder::addCTPDigit(uint32_t linkCRU, uint32_t orbit, gbtword80_t& d } } else if (linkCRU == o2::ctp::GBTLinkIDClassRec) { int32_t BCShiftCorrection = -o2::ctp::TriggerOffsetsParam::Instance().customOffset[o2::detectors::DetID::CTP]; - int32_t offset = BCShiftCorrection + o2::ctp::TriggerOffsetsParam::Instance().LM_L0 + o2::ctp::TriggerOffsetsParam::Instance().L0_L1 - 1; + int32_t offset = BCShiftCorrection + o2::ctp::TriggerOffsetsParam::Instance().LM_L0 + o2::ctp::TriggerOffsetsParam::Instance().L0_L1_classes - 1; LOG(debug) << "tcr ir ori:" << ir; if ((ir.orbit <= mTFOrbit) && ((int32_t)ir.bc < offset)) { // LOG(warning) << "Loosing tclass:" << ir; @@ -293,7 +293,12 @@ int RawDataDecoder::decodeRaw(o2::framework::InputRecord& inputs, std::vector& digitsMap, o2::pmr::vector& digits, uint32_t TFOrbit) +int RawDataDecoder::shiftInputs(std::map& digitsMap, o2::pmr::vector& digits, uint32_t TFOrbit, uint64_t trgclassmask) { // int nClasswoInp = 0; // counting classes without input which should never happen int nLM = 0; @@ -527,6 +532,7 @@ int RawDataDecoder::shiftInputs(std::map& digit int nL1 = 0; int nTwI = 0; int nTwoI = 0; + int nTwoIlost = 0; std::map digitsMapShifted; auto L0shift = o2::ctp::TriggerOffsetsParam::Instance().LM_L0; auto L1shift = L0shift + o2::ctp::TriggerOffsetsParam::Instance().L0_L1; @@ -594,11 +600,17 @@ int RawDataDecoder::shiftInputs(std::map& digit if ((d.CTPInputMask & L1MASKInputs).count()) { nL1++; } - if (d.CTPClassMask.count()) { + if ((d.CTPClassMask).to_ulong() & trgclassmask) { if (d.CTPInputMask.count()) { nTwI++; } else { - nTwoI++; + if (d.intRecord.bc == (o2::constants::lhc::LHCMaxBunches - L1shift)) { // input can be lost because latency class-l1input = 1 + nTwoIlost++; + } else { + // LOG(error) << d.intRecord << " " << d.CTPClassMask << " " << d.CTPInputMask; + // std::cout << "ERROR:" << std::hex << d.CTPClassMask << " " << d.CTPInputMask << std::dec << std::endl; + nTwoI++; + } } } digits.push_back(dig.second); @@ -606,6 +618,9 @@ int RawDataDecoder::shiftInputs(std::map& digit if (nTwoI) { // Trigger class wo Input LOG(error) << "LM:" << nLM << " L0:" << nL0 << " L1:" << nL1 << " TwI:" << nTwI << " Trigger classes wo input:" << nTwoI; } + if (nTwoIlost) { + LOG(warn) << " Trigger classes wo input from diff latency 1:" << nTwoIlost; + } return 0; } // diff --git a/Detectors/CTP/workflow/include/CTPWorkflow/EntropyDecoderSpec.h b/Detectors/CTP/workflow/include/CTPWorkflow/EntropyDecoderSpec.h index 4596fe12cb31d..eee7abb08d16c 100644 --- a/Detectors/CTP/workflow/include/CTPWorkflow/EntropyDecoderSpec.h +++ b/Detectors/CTP/workflow/include/CTPWorkflow/EntropyDecoderSpec.h @@ -34,6 +34,7 @@ class EntropyDecoderSpec : public o2::framework::Task void init(o2::framework::InitContext& ic) final; void endOfStream(o2::framework::EndOfStreamContext& ec) final; void finaliseCCDB(o2::framework::ConcreteDataMatcher& matcher, void* obj) final; + void updateTimeDependentParams(framework::ProcessingContext& pc); private: o2::ctp::CTFCoder mCTFCoder; diff --git a/Detectors/CTP/workflow/include/CTPWorkflow/RawDecoderSpec.h b/Detectors/CTP/workflow/include/CTPWorkflow/RawDecoderSpec.h index 607491b5cb48a..a5a1a75a0b594 100644 --- a/Detectors/CTP/workflow/include/CTPWorkflow/RawDecoderSpec.h +++ b/Detectors/CTP/workflow/include/CTPWorkflow/RawDecoderSpec.h @@ -16,6 +16,7 @@ #include #include "Framework/DataProcessorSpec.h" #include "Framework/Task.h" +#include "Framework/WorkflowSpec.h" #include "DataFormatsCTP/Digits.h" #include "DataFormatsCTP/LumiInfo.h" #include "CTPReconstruction/RawDataDecoder.h" @@ -50,6 +51,7 @@ class RawDecoderSpec : public framework::Task /// Input RawData: {"ROUT", "RAWDATA", 0, Lifetime::Timeframe} /// Output HW errors: {"CTP", "RAWHWERRORS", 0, Lifetime::Timeframe} -later void run(framework::ProcessingContext& ctx) final; + void updateTimeDependentParams(framework::ProcessingContext& pc); protected: private: @@ -68,6 +70,7 @@ class RawDecoderSpec : public framework::Task uint32_t mNTFToIntegrate = 1; uint32_t mNHBIntegratedT = 0; uint32_t mNHBIntegratedV = 0; + bool mDecodeinputs = 0; std::deque mHistoryT; std::deque mHistoryV; RawDataDecoder mDecoder; diff --git a/Detectors/CTP/workflow/src/EntropyDecoderSpec.cxx b/Detectors/CTP/workflow/src/EntropyDecoderSpec.cxx index 8f3da5f439f80..8c2f5d05aa031 100644 --- a/Detectors/CTP/workflow/src/EntropyDecoderSpec.cxx +++ b/Detectors/CTP/workflow/src/EntropyDecoderSpec.cxx @@ -55,9 +55,8 @@ void EntropyDecoderSpec::run(ProcessingContext& pc) mTimer.Start(false); o2::ctf::CTFIOSize iosize; - mCTFCoder.updateTimeDependentParams(pc, true); + updateTimeDependentParams(pc); auto buff = pc.inputs().get>("ctf_CTP"); - auto& digits = pc.outputs().make>(OutputRef{"digits"}); auto& lumi = pc.outputs().make(OutputRef{"CTPLumi"}); @@ -76,6 +75,20 @@ void EntropyDecoderSpec::endOfStream(EndOfStreamContext& ec) LOGF(info, "CTP Entropy Decoding total timing: Cpu: %.3e Real: %.3e s in %d slots", mTimer.CpuTime(), mTimer.RealTime(), mTimer.Counter() - 1); } +void EntropyDecoderSpec::updateTimeDependentParams(framework::ProcessingContext& pc) +{ + mCTFCoder.updateTimeDependentParams(pc, true); + if (pc.services().get().globalRunNumberChanged) { + const auto ctpcfg = pc.inputs().get("ctpconfig"); + if (mCTFCoder.getDecodeInps()) { + const auto ctpcfg = pc.inputs().get("ctpconfig"); + if (ctpcfg != nullptr) { + mCTFCoder.setCTPConfig(*ctpcfg); + LOG(info) << "ctpconfig for run done:" << mCTFCoder.getCTPConfig().getRunNumber(); + } + } + } +} DataProcessorSpec getEntropyDecoderSpec(int verbosity, unsigned int sspec) { @@ -88,7 +101,7 @@ DataProcessorSpec getEntropyDecoderSpec(int verbosity, unsigned int sspec) inputs.emplace_back("ctf_CTP", "CTP", "CTFDATA", sspec, Lifetime::Timeframe); inputs.emplace_back("ctfdict_CTP", "CTP", "CTFDICT", 0, Lifetime::Condition, ccdbParamSpec("CTP/Calib/CTFDictionaryTree")); inputs.emplace_back("trigoffset", "CTP", "Trig_Offset", 0, Lifetime::Condition, ccdbParamSpec("CTP/Config/TriggerOffsets")); - + inputs.emplace_back("ctpconfig", "CTP", "CTPCONFIG", 0, Lifetime::Condition, ccdbParamSpec("CTP/Config/Config", 1)); return DataProcessorSpec{ "ctp-entropy-decoder", inputs, diff --git a/Detectors/CTP/workflow/src/RawDecoderSpec.cxx b/Detectors/CTP/workflow/src/RawDecoderSpec.cxx index 415dbe2a1ffe3..81a927b3caee1 100644 --- a/Detectors/CTP/workflow/src/RawDecoderSpec.cxx +++ b/Detectors/CTP/workflow/src/RawDecoderSpec.cxx @@ -13,20 +13,21 @@ #include #include "Framework/InputRecordWalker.h" #include "Framework/DataRefUtils.h" -#include "Framework/WorkflowSpec.h" #include "Framework/ConfigParamRegistry.h" #include "DetectorsRaw/RDHUtils.h" #include "CTPWorkflow/RawDecoderSpec.h" #include "CommonUtils/VerbosityConfig.h" #include "Framework/InputRecord.h" #include "DataFormatsCTP/TriggerOffsetsParam.h" +#include "Framework/CCDBParamSpec.h" +#include "DataFormatsCTP/Configuration.h" using namespace o2::ctp::reco_workflow; void RawDecoderSpec::init(framework::InitContext& ctx) { - bool decodeinps = ctx.options().get("ctpinputs-decoding"); - mDecoder.setDecodeInps(decodeinps); + mDecodeinputs = ctx.options().get("ctpinputs-decoding"); + mDecoder.setDecodeInps(mDecodeinputs); mNTFToIntegrate = ctx.options().get("ntf-to-average"); mVerbose = ctx.options().get("use-verbose-mode"); int maxerrors = ctx.options().get("print-errors-num"); @@ -42,7 +43,7 @@ void RawDecoderSpec::init(framework::InitContext& ctx) mOutputLumiInfo.inp2 = inp2; mMaxInputSize = ctx.options().get("max-input-size"); mMaxInputSizeFatal = ctx.options().get("max-input-size-fatal"); - LOG(info) << "CTP reco init done. Inputs decoding here:" << decodeinps << " DoLumi:" << mDoLumi << " DoDigits:" << mDoDigits << " NTF:" << mNTFToIntegrate << " Lumi inputs:" << lumiinp1 << ":" << inp1 << " " << lumiinp2 << ":" << inp2 << " Max errors:" << maxerrors << " Max input size:" << mMaxInputSize << " MaxInputSizeFatal:" << mMaxInputSizeFatal; + LOG(info) << "CTP reco init done. Inputs decoding here:" << mDecodeinputs << " DoLumi:" << mDoLumi << " DoDigits:" << mDoDigits << " NTF:" << mNTFToIntegrate << " Lumi inputs:" << lumiinp1 << ":" << inp1 << " " << lumiinp2 << ":" << inp2 << " Max errors:" << maxerrors << " Max input size:" << mMaxInputSize << " MaxInputSizeFatal:" << mMaxInputSizeFatal; // mOutputLumiInfo.printInputs(); } void RawDecoderSpec::endOfStream(framework::EndOfStreamContext& ec) @@ -73,6 +74,7 @@ void RawDecoderSpec::endOfStream(framework::EndOfStreamContext& ec) } void RawDecoderSpec::run(framework::ProcessingContext& ctx) { + updateTimeDependentParams(ctx); mOutputDigits.clear(); std::map digits; using InputSpec = o2::framework::InputSpec; @@ -176,6 +178,7 @@ void RawDecoderSpec::run(framework::ProcessingContext& ctx) mOutputLumiInfo.orbit = lumiPointsHBF1[0].orbit; } mOutputLumiInfo.counts = mCountsT; + mOutputLumiInfo.countsFV0 = mCountsV; mOutputLumiInfo.nHBFCounted = mNHBIntegratedT; mOutputLumiInfo.nHBFCountedFV0 = mNHBIntegratedV; @@ -199,6 +202,8 @@ o2::framework::DataProcessorSpec o2::ctp::reco_workflow::getRawDecoderSpec(bool std::vector outputs; if (digits) { + inputs.emplace_back("ctpconfig", "CTP", "CTPCONFIG", 0, o2::framework::Lifetime::Condition, o2::framework::ccdbParamSpec("CTP/Config/Config", 1)); + inputs.emplace_back("trigoffset", "CTP", "Trig_Offset", 0, o2::framework::Lifetime::Condition, o2::framework::ccdbParamSpec("CTP/Config/TriggerOffsets")); outputs.emplace_back("CTP", "DIGITS", 0, o2::framework::Lifetime::Timeframe); } if (lumi) { @@ -219,3 +224,18 @@ o2::framework::DataProcessorSpec o2::ctp::reco_workflow::getRawDecoderSpec(bool {"max-input-size-fatal", o2::framework::VariantType::Bool, false, {"If true issue fatal error otherwise error on;y"}}, {"ctpinputs-decoding", o2::framework::VariantType::Bool, false, {"Inputs alignment: true - raw decoder - has to be compatible with CTF decoder: allowed options: 10,01,00"}}}}; } +void RawDecoderSpec::updateTimeDependentParams(framework::ProcessingContext& pc) +{ + if (pc.services().get().globalRunNumberChanged) { + pc.inputs().get("trigoffset"); + const auto& trigOffsParam = o2::ctp::TriggerOffsetsParam::Instance(); + LOG(info) << "updateing TroggerOffsetsParam: inputs L0_L1:" << trigOffsParam.L0_L1 << " classes L0_L1:" << trigOffsParam.L0_L1_classes; + if (mDecodeinputs) { + const auto ctpcfg = pc.inputs().get("ctpconfig"); + if (ctpcfg != nullptr) { + mDecoder.setCTPConfig(*ctpcfg); + LOG(info) << "ctpconfig for run done:" << mDecoder.getCTPConfig().getRunNumber(); + } + } + } +} diff --git a/Detectors/GLOQC/src/MatchITSTPCQC.cxx b/Detectors/GLOQC/src/MatchITSTPCQC.cxx index 56462344850d6..f0345175b9a59 100644 --- a/Detectors/GLOQC/src/MatchITSTPCQC.cxx +++ b/Detectors/GLOQC/src/MatchITSTPCQC.cxx @@ -344,10 +344,10 @@ bool MatchITSTPCQC::init() mChi2Refit = new TH1F("mChi2Refit", "Chi2 of refit; chi2", 200, 0, 300); mChi2Refit->SetOption("logy"); mChi2Refit->GetYaxis()->SetTitleOffset(1.4); - mDCAr = new TH1F("mDCAr", "DCA of TPC tracks; DCAr", 200, -100, 100); - mDCArVsPtNum = new TH2F("mDCArVsPtNum", "DCA of TPC tracks Vs Pt Num; #it{p}_{T} [GeV/c]; DCAr", 100, 0, 20., 200, -30, 30); + mDCAr = new TH1F("mDCAr", "DCA of TPC tracks; DCAr", 100, -mDCATPCCutY, mDCATPCCutY); + mDCArVsPtNum = new TH2F("mDCArVsPtNum", "DCA of TPC tracks Vs Pt Num; #it{p}_{T} [GeV/c]; DCAr", 100, 0, 20., 100, -mDCATPCCutY, mDCATPCCutY); mDCArVsPtNum->Sumw2(); - mDCArVsPtDen = new TH2F("mDCArVsPtDen", "DCA of TPC tracks Vs Pt Den; #it{p}_{T} [GeV/c]; DCAr", 100, 0, 20., 200, -30, 30); + mDCArVsPtDen = new TH2F("mDCArVsPtDen", "DCA of TPC tracks Vs Pt Den; #it{p}_{T} [GeV/c]; DCAr", 100, 0, 20., 100, -mDCATPCCutY, mDCATPCCutY); mDCArVsPtDen->Sumw2(); mFractionITSTPCmatchDCArVsPt = new TEfficiency("mFractionITSTPCmatchDCArVsPt", "Fraction of ITSTPC matched tracks wrt TPC vs DCAr; #it{p}_{T} [GeV#it{c}]; DCAr; Eff", 100, 0, 20., 200, -30, 30); diff --git a/Detectors/ITSMFT/ITS/postprocessing/studies/CMakeLists.txt b/Detectors/ITSMFT/ITS/postprocessing/studies/CMakeLists.txt index 361ab4db4fb8e..9794b69631d57 100644 --- a/Detectors/ITSMFT/ITS/postprocessing/studies/CMakeLists.txt +++ b/Detectors/ITSMFT/ITS/postprocessing/studies/CMakeLists.txt @@ -30,4 +30,6 @@ o2_target_root_dictionary(ITSPostprocessing HEADERS include/ITSStudies/ITSStudiesConfigParam.h include/ITSStudies/TrackCuts.h include/ITSStudies/TrackMethods.h -LINKDEF src/ITSStudiesLinkDef.h) \ No newline at end of file +LINKDEF src/ITSStudiesLinkDef.h) + +add_subdirectory(macros) diff --git a/Detectors/ITSMFT/ITS/postprocessing/studies/include/ITSStudies/TrackExtension.h b/Detectors/ITSMFT/ITS/postprocessing/studies/include/ITSStudies/TrackExtension.h index 2567000746559..fd5b93b0f9509 100644 --- a/Detectors/ITSMFT/ITS/postprocessing/studies/include/ITSStudies/TrackExtension.h +++ b/Detectors/ITSMFT/ITS/postprocessing/studies/include/ITSStudies/TrackExtension.h @@ -24,7 +24,7 @@ class MCKinematicsReader; namespace its::study { using mask_t = o2::dataformats::GlobalTrackID::mask_t; -o2::framework::DataProcessorSpec getTrackExtensionStudy(mask_t srcTracksMask, mask_t srcClustersMask, bool useMC, std::shared_ptr kineReader); +o2::framework::DataProcessorSpec getTrackExtensionStudy(mask_t srcTracksMask, mask_t srcClustersMask, std::shared_ptr kineReader); } // namespace its::study } // namespace o2 diff --git a/Detectors/ITSMFT/ITS/postprocessing/studies/macros/CMakeLists.txt b/Detectors/ITSMFT/ITS/postprocessing/studies/macros/CMakeLists.txt new file mode 100644 index 0000000000000..2d78e4077ec53 --- /dev/null +++ b/Detectors/ITSMFT/ITS/postprocessing/studies/macros/CMakeLists.txt @@ -0,0 +1,16 @@ +# Copyright 2019-2020 CERN and copyright holders of ALICE O2. +# See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +# All rights not expressly granted are reserved. +# +# This software is distributed under the terms of the GNU General Public +# License v3 (GPL Version 3), copied verbatim in the file "COPYING". +# +# In applying this license CERN does not waive the privileges and immunities +# granted to it by virtue of its status as an Intergovernmental Organization +# or submit itself to any jurisdiction. + +# o2_add_test_root_macro( +# PostTrackExtension.C +# PUBLIC_LINK_LIBRARIES ROOT::Hist ROOT::RIO ROOT::Core ROOT::Gpad +# LABELS its-study +# COMPILE_ONLY) diff --git a/Detectors/ITSMFT/ITS/postprocessing/studies/macros/PostTrackExtension.notest b/Detectors/ITSMFT/ITS/postprocessing/studies/macros/PostTrackExtension.notest new file mode 100644 index 0000000000000..4a7c9c4159a4b --- /dev/null +++ b/Detectors/ITSMFT/ITS/postprocessing/studies/macros/PostTrackExtension.notest @@ -0,0 +1,629 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#if !defined(__CLING__) || defined(__ROOTCLING__) +#include "TStyle.h" +#include "TFile.h" +#include "TError.h" +#include "TColor.h" +#include "TCanvas.h" +#include "TH2D.h" +#include "TF1.h" +#include "TEfficiency.h" +#include "TMarker.h" +#include "TLegend.h" +#include "TTree.h" +#include "TLatex.h" + +#include +#include +#include +#endif + +static constexpr std::array bitPatternsBefore{15, 30, 31, 60, 62, 63, 120, 124, 126}; +static constexpr std::array bitPatternsAfter{31, 47, 61, 62, 63, 79, 94, 95, 111, 121, 122, 123, 124, 125, 126, 127}; +inline bool bitsCleared(uint8_t b, uint8_t a) { return (a == b) ? true : (b & ~(a & b)) != 0; } +static constexpr std::array patternColors = { + kRed, // Red + kBlue, // Blue + kGreen, // Green + kMagenta, // Magenta + kCyan, // Cyan + kOrange, // Orange + kViolet, // Violet + kYellow, // Yellow + kPink, // Pink + kAzure, // Azure + kSpring, // Spring Green + kTeal, // Teal + kBlack, // Black + kGray, // Gray + kOrange + 7, // Light Orange + kBlue - 9 // Light Blue +}; + +// Marker styles +static constexpr std::array patternMarkers = { + 20, // Full circle + 21, // Full square + 22, // Full triangle up + 23, // Full triangle down + 24, // Open circle + 25, // Open square + 26, // Open triangle up + 27, // Open cross + 28, // Star + 29, // Plus sign + 30, // Open diamond + 31, // Full diamond + 32, // Cross + 33, // Circle with cross + 34, // X sign + 35 // Double open cross +}; + +enum Labels : unsigned int { + eAll = 0, + eGood, + eFake, + eFakeBefore, + eFakeAfter, + eFakeMix, + eTopGood, + eBotGood, + eMixGood, + eTopFake, + eBotFake, + eMixFake, + eN, +}; +static const std::array names{ + "ALL #frac{ext trks}{all trks}", + "GOOD #frac{good ext trks}{all ext trks}", + "FAKE #frac{fake trks}{all ext trks}", + "FAKE BF #frac{fake bf trks}{fake ext trks}", + "FAKE AF #frac{fake af trks}{fake ext trks}", + "FAKE MIX #frac{fake mix trks}{fake ext trks}", + // Good Top/Bot/Mix + "TOP #frac{good top ext trks}{good ext trks}", + "BOT #frac{good bot ext trks}{good ext trks}", + "MIX #frac{good mix ext trks}{good ext trks}", + // Fake Top/Bot/Mix + "TOP #frac{fake top ext trks}{fake ext trks}", + "BOT #frac{fake bot ext trks}{fake ext trks}", + "MIX #frac{fake mix ext trks}{fake ext trks}", +}; +static const std::array colors{kBlack, kGreen, kRed, kCyan, kYellow, kAzure, + // Good Top/Bot/Mix + kBlue, kOrange, kPink, + // Fake Top/Bot/Mix + kBlue, kOrange, kPink}; +static const std::array markers{20, 21, 22, 23, 27, 28, + // Good Top/Bot/Mix + 29, 33, 39, + // Fake Top/Bot/Mix + 29, 33, 39}; +static const char* const texPtX = "#it{p}_{T} (GeV/#it{c})"; +static const char* const texEff = "Efficiency"; +static const char* const texPtRes = "#sigma(#Delta#it{p}_{T}/#it{p}_{T})"; +static const char* const texDCAxyRes = "#sigma(DCA_{#it{xy}}) (#mum)"; +static const char* const texDCAzRes = "#sigma(DCA_{#it{z}}) (#mum)"; +static const char* const fitOpt{"QWMER"}; + +void setStyle(); +TEfficiency* makeEff(TFile*, const char* num, const char* den); + +template +void style(T* t, Labels lab, TLegend* leg = nullptr) +{ + t->SetMarkerStyle(markers[lab]); + t->SetMarkerColor(colors[lab]); + t->SetLineColor(colors[lab]); + if (leg) { + leg->AddEntry(t, names[lab]); + } +} + +template +void stylePattern(T* t, int i, TLegend* leg = nullptr, const char* name = nullptr) +{ + t->SetMarkerStyle(patternMarkers[i]); + t->SetMarkerColor(patternColors[i]); + t->SetLineColor(patternColors[i]); + if (leg) { + leg->AddEntry(t, name); + } +} + +void PostTrackExtension(const char* fileName = "TrackExtensionStudy.root") +{ + setStyle(); + + std::unique_ptr fIn{TFile::Open(fileName, "READ")}; + if (!fIn || fIn->IsZombie()) { + Error("", "Cannot open file %s", fileName); + return; + } + + { // Purity & Fake-Rate + auto c = new TCanvas("cPFR", "", 800, 600); + auto h = c->DrawFrame(0.05, 0.0, 10., 1.05); + h->GetXaxis()->SetTitle(texPtX); + h->GetYaxis()->SetTitle(texEff); + auto leg = new TLegend(0.35, 0.35, 0.7, 0.7); + auto eff = fIn->Get("eExtension"); + style(eff, eAll, leg); + eff->Draw("same"); + auto effPurity = fIn->Get("eExtensionPurity"); + style(effPurity, eGood, leg); + effPurity->Draw("same"); + auto effFake = fIn->Get("eExtensionFake"); + style(effFake, eFake, leg); + effFake->Draw("same"); + leg->Draw(); + gPad->SetLogx(); + gPad->SetGrid(); + c->SaveAs("trkExt_purity_fake.pdf"); + } + + { // FAKE-Rate composition + auto c = new TCanvas("cFR", "", 800, 600); + auto h = c->DrawFrame(0.05, 0.0, 10., 1.05); + h->GetXaxis()->SetTitle(texPtX); + h->GetYaxis()->SetTitle(texEff); + auto leg = new TLegend(0.35, 0.35, 0.7, 0.7); + auto effFake = fIn->Get("eExtensionFake"); + style(effFake, eFake, leg); + effFake->Draw("same"); + auto effFakeBf = fIn->Get("eExtensionFakeBefore"); + style(effFakeBf, eFakeBefore, leg); + effFakeBf->Draw("same"); + auto effFakeAf = fIn->Get("eExtensionFakeAfter"); + style(effFakeAf, eFakeAfter, leg); + effFakeAf->Draw("same"); + auto effFakeMi = fIn->Get("eExtensionFakeMix"); + style(effFakeMi, eFakeMix, leg); + effFakeMi->Draw("same"); + leg->Draw(); + gPad->SetLogx(); + gPad->SetGrid(); + c->SaveAs("trkExt_fake.pdf"); + } + + { // GOOD Top/Bot/Mix Purity composition + auto c = new TCanvas("cGC", "", 800, 600); + auto h = c->DrawFrame(0.05, 0.0, 10., 1.05); + h->GetXaxis()->SetTitle(texPtX); + h->GetYaxis()->SetTitle(texEff); + auto leg = new TLegend(0.35, 0.35, 0.7, 0.7); + auto effTop = makeEff(fIn.get(), "eExtensionTopPurity", "eExtensionPurity"); + style(effTop, eTopGood, leg); + effTop->Draw("same"); + auto effBot = makeEff(fIn.get(), "eExtensionBotPurity", "eExtensionPurity"); + style(effBot, eBotGood, leg); + effBot->Draw("same"); + auto effMix = makeEff(fIn.get(), "eExtensionMixPurity", "eExtensionPurity"); + style(effMix, eMixGood, leg); + effMix->Draw("same"); + leg->Draw(); + gPad->SetLogx(); + gPad->SetGrid(); + c->SaveAs("trkExt_good_comp.pdf"); + } + + { // FAKE Top/Bot/Mix composition + auto c = new TCanvas("cFC", "", 800, 600); + auto h = c->DrawFrame(0.05, 0.0, 10., 1.05); + h->GetXaxis()->SetTitle(texPtX); + h->GetYaxis()->SetTitle(texEff); + auto leg = new TLegend(0.35, 0.35, 0.7, 0.7); + auto effTop = fIn->Get("eExtensionTopFake"); + style(effTop, eTopFake, leg); + effTop->Draw("same"); + auto effBot = fIn->Get("eExtensionBotFake"); + style(effBot, eBotFake, leg); + effBot->Draw("same"); + auto effMix = fIn->Get("eExtensionMixFake"); + style(effMix, eMixFake, leg); + effMix->Draw("same"); + leg->Draw(); + gPad->SetLogx(); + gPad->SetGrid(); + c->SaveAs("trkExt_fake_comp.pdf"); + } + + { // Good Patterns + auto c = new TCanvas("cPatGood", "", 3 * 800, 3 * 600); + c->Divide(3, 3); + for (int i{0}; i < (int)bitPatternsBefore.size(); ++i) { + auto p = c->cd(i + 1); + auto h = p->DrawFrame(0.05, 0.0, 10., 1.05); + h->GetXaxis()->SetTitle(texPtX); + h->GetYaxis()->SetTitle(texEff); + auto leg = new TLegend(0.35, 0.60, 0.7, 0.88); + leg->SetNColumns(4); + leg->SetHeader(std::format("BEFORE={:07b} GOOD Pattern AFTER/BEFORE", bitPatternsBefore[i]).c_str()); + for (int j{0}; j < (int)bitPatternsAfter.size(); ++j) { + if (bitsCleared(bitPatternsBefore[i], bitPatternsAfter[j])) { + continue; + } + auto eff = fIn->Get(std::format("eExtensionPatternGood_{:07b}_{:07b}", bitPatternsBefore[i], bitPatternsAfter[j]).c_str()); + stylePattern(eff, j, leg, std::format("{:07b}", bitPatternsAfter[j]).c_str()); + eff->Draw("same"); + } + leg->Draw(); + p->SetLogx(); + p->SetGrid(); + } + c->SaveAs("trkExt_good_pattern_comp.pdf"); + } + + { // Fake Patterns + auto c = new TCanvas("cPatFake", "", 3 * 800, 3 * 600); + c->Divide(3, 3); + for (int i{0}; i < (int)bitPatternsBefore.size(); ++i) { + auto p = c->cd(i + 1); + auto h = p->DrawFrame(0.05, 0.0, 10., 1.05); + h->GetXaxis()->SetTitle(texPtX); + h->GetYaxis()->SetTitle(texEff); + auto leg = new TLegend(0.35, 0.60, 0.7, 0.88); + leg->SetNColumns(4); + leg->SetHeader(std::format("BEFORE={:07b} FAKE Pattern AFTER/BEFORE", bitPatternsBefore[i]).c_str()); + for (int j{0}; j < (int)bitPatternsAfter.size(); ++j) { + if (bitsCleared(bitPatternsBefore[i], bitPatternsAfter[j])) { + continue; + } + auto eff = fIn->Get(std::format("eExtensionPatternFake_{:07b}_{:07b}", bitPatternsBefore[i], bitPatternsAfter[j]).c_str()); + stylePattern(eff, j, leg, std::format("{:07b}", bitPatternsAfter[j]).c_str()); + eff->Draw("same"); + } + leg->Draw(); + p->SetLogx(); + p->SetGrid(); + } + c->SaveAs("trkExt_fake_pattern_comp.pdf"); + } + + { // DCA + auto fGaus = new TF1("fGaus", "gaus", -200., 200.); + auto dcaXYVsPtNo = fIn->Get("hDCAxyVsPtResNormal"); + auto dcaXYVsPtYes = fIn->Get("hDCAxyVsPtResExtended"); + auto dcazVsPtNo = fIn->Get("hDCAzVsPtResNormal"); + auto dcazVsPtYes = fIn->Get("hDCAzVsPtResExtended"); + auto bins = dcazVsPtNo->GetXaxis()->GetXbins(); + auto dcaXYResNo = new TH1F("hDcaxyResNo", "NORMAL;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", bins->GetSize() - 1, bins->GetArray()); + auto dcaXYResYes = new TH1F("hDcaxyResYes", "EXTENDED;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", bins->GetSize() - 1, bins->GetArray()); + auto dcaZResNo = new TH1F("hDcazResNo", "NORMAL;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", bins->GetSize() - 1, bins->GetArray()); + auto dcaZResYes = new TH1F("hDcazResYes", "EXTENDED;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", bins->GetSize() - 1, bins->GetArray()); + TH1* proj; + for (int iPt{1}; iPt <= bins->GetSize(); ++iPt) { + auto ptMin = dcaXYResNo->GetXaxis()->GetBinLowEdge(iPt); + if (ptMin < 0.1) { + continue; + } + float minFit = (ptMin < 1.) ? -200. : -75.; + float maxFit = (ptMin < 1.) ? 200. : 75.; + + proj = dcaXYVsPtNo->ProjectionY(Form("hProjDCAxy_no_%d", iPt), iPt, iPt); + proj->Fit("fGaus", fitOpt, "", minFit, maxFit); + dcaXYResNo->SetBinContent(iPt, fGaus->GetParameter(2)); + dcaXYResNo->SetBinError(iPt, fGaus->GetParError(2)); + + proj = dcaXYVsPtYes->ProjectionY(Form("hProjDCAxy_yes_%d", iPt), iPt, iPt); + proj->Fit("fGaus", fitOpt, "", minFit, maxFit); + dcaXYResYes->SetBinContent(iPt, fGaus->GetParameter(2)); + dcaXYResYes->SetBinError(iPt, fGaus->GetParError(2)); + + proj = dcazVsPtNo->ProjectionY(Form("hProjDCAz_no_%d", iPt), iPt, iPt); + proj->Fit("fGaus", fitOpt, "", minFit, maxFit); + dcaZResNo->SetBinContent(iPt, fGaus->GetParameter(2)); + dcaZResNo->SetBinError(iPt, fGaus->GetParError(2)); + + proj = dcazVsPtYes->ProjectionY(Form("hProjDCAz_yes_%d", iPt), iPt, iPt); + proj->Fit("fGaus", fitOpt, "", minFit, maxFit); + dcaZResYes->SetBinContent(iPt, fGaus->GetParameter(2)); + dcaZResYes->SetBinError(iPt, fGaus->GetParError(2)); + } + + dcaXYResNo->SetLineColor(kRed); + dcaXYResNo->SetMarkerColor(kRed); + dcaXYResYes->SetLineColor(kBlue); + dcaXYResYes->SetMarkerColor(kBlue); + dcaZResNo->SetLineColor(kRed); + dcaZResNo->SetMarkerColor(kRed); + dcaZResYes->SetLineColor(kBlue); + dcaZResYes->SetMarkerColor(kBlue); + + auto c = new TCanvas("cDCA", "", 2 * 800, 600); + c->Divide(2, 1); + c->cd(1); + auto h = gPad->DrawFrame(0.1, 1, 10., 500); + h->GetXaxis()->SetTitle(texPtX); + h->GetYaxis()->SetTitle(texDCAxyRes); + dcaXYResNo->Draw("SAME"); + dcaXYResYes->Draw("SAME"); + gPad->SetLogx(); + gPad->SetLogy(); + gPad->SetGrid(); + auto leg = new TLegend(0.20, 0.20, 0.40, 0.40); + leg->AddEntry(dcaXYResNo, "Normal"); + leg->AddEntry(dcaXYResYes, "Extended"); + leg->Draw(); + + c->cd(2); + h = gPad->DrawFrame(0.1, 1, 10., 500); + h->GetXaxis()->SetTitle(texPtX); + h->GetYaxis()->SetTitle(texDCAzRes); + dcaZResNo->Draw("SAME"); + dcaZResYes->Draw("SAME"); + gPad->SetLogx(); + gPad->SetLogy(); + gPad->SetGrid(); + + c->SaveAs("trkExt_dca.pdf"); + } + + return; + { // Kinematic variables + auto t = fIn->Get("tree"); + auto c = new TCanvas("cKG", "", 800, 600); + c->Divide(3, 2); + { + auto p = c->cd(1); + p->SetGrid(); + auto h = p->DrawFrame(-.6, 0., .6, 9.); + h->GetXaxis()->SetTitle("#frac{Q^{2}}{p_{T,TRK}}-#frac{Q^{2}}{p_{T,MC}}"); + h->GetYaxis()->SetTitle("n. counts"); + t->Draw("trk.getQ2Pt()-mcTrk.getQ2Pt()>>hPtNo(100,-.6,.6)", "isGood&&!isExtended", "HIST;SAME"); + auto hNo = (TH1F*)p->GetPrimitive("hPtNo"); + hNo->Scale(1.0 / hNo->Integral("width")); + hNo->SetLineColor(kRed); + auto fitNo = new TF1("fitNo", "gaus", -0.04, 0.04); + hNo->Fit(fitNo, "QR"); + fitNo->SetLineColor(kRed); + fitNo->Draw("SAME"); + auto textNo = new TLatex(-0.55, 8.2, Form("#mu = %.3f, #sigma = %.3f", fitNo->GetParameter(1), fitNo->GetParameter(2))); + textNo->SetTextColor(kRed); + textNo->SetNDC(false); + textNo->SetTextSize(0.05); + textNo->Draw(); + + t->Draw("trk.getQ2Pt()-mcTrk.getQ2Pt()>>hPtYes(100,-.6,.6)", "isGood&&isExtended", "HIST;SAME"); + auto hYes = (TH1F*)p->GetPrimitive("hPtYes"); + hYes->Scale(1.0 / hYes->Integral("width")); + hYes->SetLineColor(kBlue); + auto fitYes = new TF1("fitYes", "gaus", -0.04, 0.04); + hYes->Fit(fitYes, "QR"); + fitYes->SetLineColor(kBlue); + fitYes->Draw("SAME"); + auto textYes = new TLatex(-0.55, 7, Form("#mu = %.4f, #sigma = %.4f", fitNo->GetParameter(1), fitNo->GetParameter(2))); + textYes->SetTextColor(kBlue); + textYes->SetNDC(false); + textYes->SetTextSize(0.05); + textYes->Draw(); + + p->Modified(); + p->Update(); + } + { + auto p = c->cd(2); + p->SetGrid(); + auto h = p->DrawFrame(-3, 0., 3, 2.); + h->GetXaxis()->SetTitle("Y_{TRK}-Y_{MC}"); + h->GetYaxis()->SetTitle("n. counts"); + t->Draw("trk.getY()-mcTrk.getY()>>hYNo(100,-3,3)", "isGood&&!isExtended", "HIST;SAME"); + auto hNo = (TH1F*)p->GetPrimitive("hYNo"); + hNo->Scale(1.0 / hNo->Integral("width")); + hNo->SetLineColor(kRed); + auto fitNo = new TF1("fitNo", "gaus", -0.5, 0.5); + hNo->Fit(fitNo, "QR"); + fitNo->SetLineColor(kRed); + fitNo->Draw("SAME"); + auto textNo = new TLatex(-2, 1.7, Form("#mu = %.3f, #sigma = %.3f", fitNo->GetParameter(1), fitNo->GetParameter(2))); + textNo->SetTextColor(kRed); + textNo->SetNDC(false); + textNo->SetTextSize(0.05); + textNo->Draw(); + + t->Draw("trk.getY()-mcTrk.getY()>>hYYes(100,-3,3)", "isGood&&isExtended", "HIST;SAME"); + auto hYes = (TH1F*)p->GetPrimitive("hYYes"); + hYes->Scale(1.0 / hYes->Integral("width")); + hYes->SetLineColor(kBlue); + auto fitYes = new TF1("fitYes", "gaus", -0.5, 0.5); + hYes->Fit(fitYes, "QR"); + fitYes->SetLineColor(kBlue); + fitYes->Draw("SAME"); + auto textYes = new TLatex(-2, 1.5, Form("#mu = %.4f, #sigma = %.4f", fitNo->GetParameter(1), fitNo->GetParameter(2))); + textYes->SetTextColor(kBlue); + textYes->SetNDC(false); + textYes->SetTextSize(0.05); + textYes->Draw(); + + p->Modified(); + p->Update(); + } + { + auto p = c->cd(3); + p->SetGrid(); + auto h = p->DrawFrame(-2, 0., 2, 4.2); + h->GetXaxis()->SetTitle("Z_{TRK}-Z_{MC}"); + h->GetYaxis()->SetTitle("n. counts"); + t->Draw("trk.getZ()-mcTrk.getZ()>>hZNo(100,-2,2)", "isGood&&!isExtended", "HIST;SAME"); + auto hNo = (TH1F*)p->GetPrimitive("hZNo"); + hNo->Scale(1.0 / hNo->Integral("width")); + hNo->SetLineColor(kRed); + auto fitNo = new TF1("fitNo", "gaus", -0.2, 0.2); + hNo->Fit(fitNo, "QR"); + fitNo->SetLineColor(kRed); + fitNo->Draw("SAME"); + auto textNo = new TLatex(-1.7, 3.8, Form("#mu = %.3f, #sigma = %.3f", fitNo->GetParameter(1), fitNo->GetParameter(2))); + textNo->SetTextColor(kRed); + textNo->SetNDC(false); + textNo->SetTextSize(0.05); + textNo->Draw(); + + t->Draw("trk.getZ()-mcTrk.getZ()>>hZYes(100,-2,2)", "isGood&&isExtended", "HIST;SAME"); + auto hYes = (TH1F*)p->GetPrimitive("hZYes"); + hYes->Scale(1.0 / hYes->Integral("width")); + hYes->SetLineColor(kBlue); + auto fitYes = new TF1("fitYes", "gaus", -0.2, 0.2); + hYes->Fit(fitYes, "QR"); + fitYes->SetLineColor(kBlue); + fitYes->Draw("SAME"); + auto textYes = new TLatex(-1.7, 3.5, Form("#mu = %.4f, #sigma = %.4f", fitNo->GetParameter(1), fitNo->GetParameter(2))); + textYes->SetTextColor(kBlue); + textYes->SetNDC(false); + textYes->SetTextSize(0.05); + textYes->Draw(); + + p->Modified(); + p->Update(); + } + { + auto p = c->cd(4); + p->SetGrid(); + auto h = p->DrawFrame(-0.02, 0., 0.02, 370.); + h->GetXaxis()->SetTitle("TGL_{TRK}-TGL_{MC}"); + h->GetYaxis()->SetTitle("n. counts"); + t->Draw("trk.getTgl()-mcTrk.getTgl()>>hTglNo(100,-0.02,0.02)", "isGood&&!isExtended", "HIST;SAME"); + auto hNo = (TH1F*)p->GetPrimitive("hTglNo"); + hNo->Scale(1.0 / hNo->Integral("width")); + hNo->SetLineColor(kRed); + auto fitNo = new TF1("fitNo", "gaus", -0.003, 0.003); + hNo->Fit(fitNo, "QR"); + fitNo->SetLineColor(kRed); + fitNo->Draw("SAME"); + auto textNo = new TLatex(-0.018, 330, Form("#mu = %.3f, #sigma = %.3f", fitNo->GetParameter(1), fitNo->GetParameter(2))); + textNo->SetTextColor(kRed); + textNo->SetNDC(false); + textNo->SetTextSize(0.05); + textNo->Draw(); + + t->Draw("trk.getTgl()-mcTrk.getTgl()>>hTglYes(100,-0.02,0.02)", "isGood&&isExtended", "HIST;SAME"); + auto hYes = (TH1F*)p->GetPrimitive("hTglYes"); + hYes->Scale(1.0 / hYes->Integral("width")); + hYes->SetLineColor(kBlue); + auto fitYes = new TF1("fitYes", "gaus", -0.003, 0.003); + hYes->Fit(fitYes, "QR"); + fitYes->SetLineColor(kBlue); + fitYes->Draw("SAME"); + auto textYes = new TLatex(-0.018, 310, Form("#mu = %.6f, #sigma = %.6f", fitNo->GetParameter(1), fitNo->GetParameter(2))); + textYes->SetTextColor(kBlue); + textYes->SetNDC(false); + textYes->SetTextSize(0.05); + textYes->Draw(); + + p->Modified(); + p->Update(); + } + { + auto p = c->cd(5); + p->SetGrid(); + auto h = p->DrawFrame(-0.08, 0., 0.08, 80.); + h->GetXaxis()->SetTitle("SNP_{TRK}-SNP_{MC}"); + h->GetYaxis()->SetTitle("n. counts"); + t->Draw("trk.getSnp()-mcTrk.getSnp()>>hSnpNo(100,-0.08,0.08)", "isGood&&!isExtended", "HIST;SAME"); + auto hNo = (TH1F*)p->GetPrimitive("hSnpNo"); + hNo->Scale(1.0 / hNo->Integral("width")); + hNo->SetLineColor(kRed); + auto fitNo = new TF1("fitNo", "gaus", -0.03, 0.03); + hNo->Fit(fitNo, "QR"); + fitNo->SetLineColor(kRed); + fitNo->Draw("SAME"); + auto textNo = new TLatex(-0.07, 72, Form("#mu = %.3f, #sigma = %.3f", fitNo->GetParameter(1), fitNo->GetParameter(2))); + textNo->SetTextColor(kRed); + textNo->SetNDC(false); + textNo->SetTextSize(0.05); + textNo->Draw(); + + t->Draw("trk.getSnp()-mcTrk.getSnp()>>hSnpYes(100,-0.08,0.08)", "isGood&&isExtended", "HIST;SAME"); + auto hYes = (TH1F*)p->GetPrimitive("hSnpYes"); + hYes->Scale(1.0 / hYes->Integral("width")); + hYes->SetLineColor(kBlue); + auto fitYes = new TF1("fitYes", "gaus", -0.03, 0.03); + hYes->Fit(fitYes, "QR"); + fitYes->SetLineColor(kBlue); + fitYes->Draw("SAME"); + auto textYes = new TLatex(-0.07, 66, Form("#mu = %.6f, #sigma = %.6f", fitNo->GetParameter(1), fitNo->GetParameter(2))); + textYes->SetTextColor(kBlue); + textYes->SetNDC(false); + textYes->SetTextSize(0.05); + textYes->Draw(); + + p->Modified(); + p->Update(); + } + { + auto p = c->cd(6); + auto legend = new TLegend(0.2, 0.2, 0.8, 0.8); + legend->SetTextSize(0.06); + legend->SetLineWidth(3); + legend->SetHeader("GOOD tracks", "C"); + auto mBlue = new TMarker(); + mBlue->SetMarkerColor(kBlue); + mBlue->SetMarkerSize(4); + legend->AddEntry(mBlue, "extended", "p"); + auto mRed = new TMarker(); + mRed->SetMarkerColor(kRed); + mRed->SetMarkerSize(4); + legend->AddEntry(mRed, "normal", "p"); + legend->SetLineColor(kRed); + legend->Draw(); + } + c->SaveAs("trkExt_kinematics.pdf"); + } +} + +void setStyle() +{ + gStyle->Reset("Plain"); + gStyle->SetOptTitle(0); + gStyle->SetOptStat(0); + gStyle->SetPalette(kRainbow); + gStyle->SetCanvasColor(10); + gStyle->SetCanvasBorderMode(0); + gStyle->SetFrameLineWidth(1); + gStyle->SetFrameFillColor(kWhite); + gStyle->SetPadColor(10); + gStyle->SetPadTickX(1); + gStyle->SetPadTickY(1); + gStyle->SetPadBottomMargin(0.15); + gStyle->SetPadLeftMargin(0.15); + gStyle->SetHistLineWidth(1); + gStyle->SetHistLineColor(kRed); + gStyle->SetFuncWidth(2); + gStyle->SetFuncColor(kGreen); + gStyle->SetLineWidth(2); + gStyle->SetLabelSize(0.045, "xyz"); + gStyle->SetLabelOffset(0.01, "y"); + gStyle->SetLabelOffset(0.01, "x"); + gStyle->SetLabelColor(kBlack, "xyz"); + gStyle->SetTitleSize(0.05, "xyz"); + gStyle->SetTitleOffset(1.25, "y"); + gStyle->SetTitleOffset(1.2, "x"); + gStyle->SetTitleFillColor(kWhite); + gStyle->SetTextSizePixels(26); + gStyle->SetTextFont(42); + gStyle->SetTickLength(0.04, "X"); + gStyle->SetTickLength(0.04, "Y"); + gStyle->SetLegendBorderSize(0); + gStyle->SetLegendFillColor(kWhite); + gStyle->SetFillColor(kWhite); + gStyle->SetLegendFont(42); +} + +TEfficiency* makeEff(TFile* fIn, const char* num, const char* den) +{ + auto h1 = fIn->Get(num)->GetPassedHistogram(); + auto h2 = fIn->Get(den)->GetPassedHistogram(); + auto e = new TEfficiency(*h1, *h2); + return e; +} diff --git a/Detectors/ITSMFT/ITS/postprocessing/studies/src/TrackExtension.cxx b/Detectors/ITSMFT/ITS/postprocessing/studies/src/TrackExtension.cxx index 364a354c700b6..465365ffa3d86 100644 --- a/Detectors/ITSMFT/ITS/postprocessing/studies/src/TrackExtension.cxx +++ b/Detectors/ITSMFT/ITS/postprocessing/studies/src/TrackExtension.cxx @@ -14,17 +14,24 @@ #include "DataFormatsITS/TrackITS.h" #include "DataFormatsITSMFT/CompCluster.h" #include "DetectorsBase/GRPGeomHelper.h" +#include "DetectorsBase/Propagator.h" #include "Framework/ConfigParamRegistry.h" #include "Framework/Task.h" #include "ITSBase/GeometryTGeo.h" #include "ITSStudies/Helpers.h" #include "ITSStudies/TrackExtension.h" +#include "SimulationDataFormat/MCEventHeader.h" #include "SimulationDataFormat/MCTrack.h" #include "Steer/MCKinematicsReader.h" +#include "ReconstructionDataFormats/Vertex.h" +#include "ReconstructionDataFormats/DCA.h" #include #include "TFile.h" +#include "TH1D.h" +#include "TH2D.h" +#include "TEfficiency.h" namespace o2::its::study { @@ -36,7 +43,9 @@ using o2::steer::MCKinematicsReader; class TrackExtensionStudy : public Task { struct ParticleInfo { - int event; + float eventX; + float eventY; + float eventZ; int pdg; float pt; float eta; @@ -60,24 +69,24 @@ class TrackExtensionStudy : public Task public: TrackExtensionStudy(std::shared_ptr dr, mask_t src, - bool useMC, std::shared_ptr kineReader, std::shared_ptr gr) : mDataRequest(dr), mTracksSrc(src), mKineReader(kineReader), mGGCCDBRequest(gr) { - if (useMC) { - LOGP(info, "Read MCKine reader with {} sources", mKineReader->getNSources()); - } + LOGP(info, "Read MCKine reader with {} sources", mKineReader->getNSources()); } ~TrackExtensionStudy() final = default; void init(InitContext& /*ic*/) final; void run(ProcessingContext& /*pc*/) final; void endOfStream(EndOfStreamContext& /*ec*/) final; + void finaliseCCDB(ConcreteDataMatcher& matcher, void* obj) final; void process(); private: static constexpr std::array mBitPatternsBefore{15, 30, 31, 60, 62, 63, 120, 124, 126}; static constexpr std::array mBitPatternsAfter{31, 47, 61, 62, 63, 79, 94, 95, 111, 121, 122, 123, 124, 125, 126, 127}; + const std::bitset<7> mTopMask{"1110000"}; + const std::bitset<7> mBotMask{"0000111"}; void updateTimeDependentParams(ProcessingContext& pc); std::string mOutFileName = "TrackExtensionStudy.root"; @@ -101,12 +110,30 @@ class TrackExtensionStudy : public Task bool mWithTree{false}; std::unique_ptr mHTrackCounts; - std::unique_ptr mHClustersCounts; std::unique_ptr mHLengthAny, mHLengthGood, mHLengthFake; std::unique_ptr mHChi2Any, mHChi2Good, mHChi2Fake; std::unique_ptr mHPtAny, mHPtGood, mHPtFake; std::unique_ptr mHExtensionAny, mHExtensionGood, mHExtensionFake; - std::unique_ptr mHExtensionPatternsAny, mHExtensionPatternsGood, mHExtensionPatternsFake; + std::unique_ptr mHExtensionPatternsAny, mHExtensionPatternsGood, mHExtensionPatternsFake, mHExtensionPatternsGoodMissed, mHExtensionPatternsGoodEmpty; + std::unique_ptr mEExtensionNum, mEExtensionDen, mEExtensionPurityNum, mEExtensionPurityDen, mEExtensionFakeNum, mEExtensionFakeDen; + std::unique_ptr mEExtensionFakeBeforeNum, mEExtensionFakeAfterNum, mEExtensionFakeMixNum; + std::unique_ptr mEExtensionTopNum, mEExtensionTopPurityNum, mEExtensionTopFakeNum; + std::unique_ptr mEExtensionBotNum, mEExtensionBotPurityNum, mEExtensionBotFakeNum; + std::unique_ptr mEExtensionMixNum, mEExtensionMixPurityNum, mEExtensionMixFakeNum; + std::array, mBitPatternsBefore.size()> mEExtensionPatternGoodNum, mEExtensionPatternFakeNum; + std::array, mBitPatternsAfter.size()>, mBitPatternsBefore.size()> mEExtensionPatternIndGoodNum, mEExtensionPatternIndFakeNum; + // DCA + std::unique_ptr mDCAxyVsPtPionsNormal, mDCAxyVsPtPionsExtended; + std::unique_ptr mDCAzVsPtPionsNormal, mDCAzVsPtPionsExtended; + + template + std::unique_ptr createHistogram(C... n, F... b) + { + auto t = std::make_unique(n..., b...); + mHistograms.push_back(static_cast(t.get())); + return std::move(t); + } + std::vector mHistograms; }; void TrackExtensionStudy::init(InitContext& ic) @@ -114,13 +141,13 @@ void TrackExtensionStudy::init(InitContext& ic) o2::base::GRPGeomHelper::instance().setRequest(mGGCCDBRequest); mWithTree = ic.options().get("with-tree"); - constexpr size_t effHistBins = 100; + constexpr size_t effHistBins = 40; constexpr float effPtCutLow = 0.01; constexpr float effPtCutHigh = 10.; auto xbins = helpers::makeLogBinning(effHistBins, effPtCutLow, effPtCutHigh); // Track Counting - mHTrackCounts = std::make_unique("hTrackCounts", "Track Stats", 10, 0, 10); + mHTrackCounts = createHistogram("hTrackCounts", "Track Stats", 10, 0, 10); mHTrackCounts->GetXaxis()->SetBinLabel(1, "Total Tracks"); mHTrackCounts->GetXaxis()->SetBinLabel(2, "Normal ANY Tracks"); mHTrackCounts->GetXaxis()->SetBinLabel(3, "Normal GOOD Tracks"); @@ -132,49 +159,87 @@ void TrackExtensionStudy::init(InitContext& ic) mHTrackCounts->GetXaxis()->SetBinLabel(9, "Extended FAKE AFTER Tracks"); mHTrackCounts->GetXaxis()->SetBinLabel(10, "Extended FAKE BEFORE&AFTER Tracks"); - // Cluster Counting - mHClustersCounts = std::make_unique("hClusterCounts", "Cluster Stats", 5, 0, 5); - mHClustersCounts->GetXaxis()->SetBinLabel(1, "Total Clusters"); - mHClustersCounts->GetXaxis()->SetBinLabel(2, "Tracking"); - mHClustersCounts->GetXaxis()->SetBinLabel(3, "Extension"); - mHClustersCounts->GetXaxis()->SetBinLabel(4, "Good Extension"); - mHClustersCounts->GetXaxis()->SetBinLabel(5, "Fake Extension"); - // Length - mHLengthAny = std::make_unique("hLengthAny", "Extended Tracks Length (ANY);NCluster;Entries", 5, 3, 8); - mHLengthGood = std::make_unique("hLengthGood", "Extended Tracks Length (GOOD);NCluster;Entries", 5, 3, 8); - mHLengthFake = std::make_unique("hLengthFake", "Extended Tracks Length (FAKE);NCluster;Entries", 5, 3, 8); + mHLengthAny = createHistogram("hLengthAny", "Extended Tracks Length (ANY);NCluster;Entries", 5, 3, 8); + mHLengthGood = createHistogram("hLengthGood", "Extended Tracks Length (GOOD);NCluster;Entries", 5, 3, 8); + mHLengthFake = createHistogram("hLengthFake", "Extended Tracks Length (FAKE);NCluster;Entries", 5, 3, 8); // Chi2 - mHChi2Any = std::make_unique("hChi2Any", "Extended Tracks Length (ANY);#chi^{2};Entries", 50, 0, 100); - mHChi2Good = std::make_unique("hChi2Good", "Extended Tracks Length (GOOD);#chi^{2};Entries", 50, 0, 100); - mHChi2Fake = std::make_unique("hChi2Fake", "Extended Tracks Length (FAKE);#chi^{2};Entries", 50, 0, 100); + mHChi2Any = createHistogram("hChi2Any", "Extended Tracks Length (ANY);#chi^{2};Entries", 50, 0, 100); + mHChi2Good = createHistogram("hChi2Good", "Extended Tracks Length (GOOD);#chi^{2};Entries", 50, 0, 100); + mHChi2Fake = createHistogram("hChi2Fake", "Extended Tracks Length (FAKE);#chi^{2};Entries", 50, 0, 100); // Pt - mHPtAny = std::make_unique("hPtAny", "Extended Tracks Length (ANY);#it{p}_{T};Entries", xbins.size(), effPtCutLow, effPtCutHigh); - mHPtGood = std::make_unique("hPtGood", "Extended Tracks Length (GOOD);#it{p}_{T};Entries", xbins.size(), effPtCutLow, effPtCutHigh); - mHPtFake = std::make_unique("hPtFake", "Extended Tracks Length (FAKE);#it{p}_{T};Entries", xbins.size(), effPtCutLow, effPtCutHigh); + mHPtAny = createHistogram("hPtAny", "Extended Tracks Length (ANY);#it{p}_{T};Entries", effHistBins, xbins.data()); + mHPtGood = createHistogram("hPtGood", "Extended Tracks Length (GOOD);#it{p}_{T};Entries", effHistBins, xbins.data()); + mHPtFake = createHistogram("hPtFake", "Extended Tracks Length (FAKE);#it{p}_{T};Entries", effHistBins, xbins.data()); // Length - mHExtensionAny = std::make_unique("hExtensionAny", "Extended Tracks Length (ANY);Extended Layer;Entries", 7, 0, 7); - mHExtensionGood = std::make_unique("hExtensionGood", "Extended Tracks Length (GOOD);Extended Layer;Entries", 7, 0, 7); - mHExtensionFake = std::make_unique("hExtensionFake", "Extended Tracks Length (FAKE);Extended Layer;Entries", 7, 0, 7); + mHExtensionAny = createHistogram("hExtensionAny", "Extended Tracks Length (ANY);Extended Layer;Entries", 7, 0, 7); + mHExtensionGood = createHistogram("hExtensionGood", "Extended Tracks Length (GOOD);Extended Layer;Entries", 7, 0, 7); + mHExtensionFake = createHistogram("hExtensionFake", "Extended Tracks Length (FAKE);Extended Layer;Entries", 7, 0, 7); // Patterns - auto makePatternAxisLabels = [&](TH1* h) { + auto makePatternAxisLabels = [&](TH1* h, bool xBefore = true) { for (int i{1}; i <= h->GetXaxis()->GetNbins(); ++i) { - h->GetXaxis()->SetBinLabel(i, fmt::format("{:07b}", mBitPatternsBefore[i - 1]).c_str()); + if (xBefore) { + h->GetXaxis()->SetBinLabel(i, fmt::format("{:07b}", mBitPatternsBefore[i - 1]).c_str()); + } else { + h->GetXaxis()->SetBinLabel(i, fmt::format("{:07b}", mBitPatternsAfter[i - 1]).c_str()); + } } for (int i{1}; i <= h->GetYaxis()->GetNbins(); ++i) { h->GetYaxis()->SetBinLabel(i, fmt::format("{:07b}", mBitPatternsAfter[i - 1]).c_str()); } }; - mHExtensionPatternsAny = std::make_unique("hExtensionPatternsAny", "Extended Tracks Pattern (ANY);Before;After;Entries", mBitPatternsBefore.size(), 0, mBitPatternsBefore.size(), mBitPatternsAfter.size(), 0, mBitPatternsAfter.size()); + mHExtensionPatternsAny = createHistogram("hExtensionPatternsAny", "Extended Tracks Pattern (ANY);Before;After;Entries", mBitPatternsBefore.size(), 0, mBitPatternsBefore.size(), mBitPatternsAfter.size(), 0, mBitPatternsAfter.size()); makePatternAxisLabels(mHExtensionPatternsAny.get()); - mHExtensionPatternsGood = std::make_unique("hExtensionPatternsGood", "Extended Tracks Pattern (GOOD);Before;After;Entries", mBitPatternsBefore.size(), 0, mBitPatternsBefore.size(), mBitPatternsAfter.size(), 0, mBitPatternsAfter.size()); + mHExtensionPatternsGood = createHistogram("hExtensionPatternsGood", "Extended Tracks Pattern (GOOD);Before;After;Entries", mBitPatternsBefore.size(), 0, mBitPatternsBefore.size(), mBitPatternsAfter.size(), 0, mBitPatternsAfter.size()); makePatternAxisLabels(mHExtensionPatternsGood.get()); - mHExtensionPatternsFake = std::make_unique("hExtensionPatternsFake", "Extended Tracks Pattern (FAKE);Before;After;Entries", mBitPatternsBefore.size(), 0, mBitPatternsBefore.size(), mBitPatternsAfter.size(), 0, mBitPatternsAfter.size()); + mHExtensionPatternsFake = createHistogram("hExtensionPatternsFake", "Extended Tracks Pattern (FAKE);Before;After;Entries", mBitPatternsBefore.size(), 0, mBitPatternsBefore.size(), mBitPatternsAfter.size(), 0, mBitPatternsAfter.size()); makePatternAxisLabels(mHExtensionPatternsFake.get()); + mHExtensionPatternsGoodMissed = createHistogram("hExtensionPatternsGoodMissed", "Extended Tracks Pattern (GOOD) Missed Clusters;After;Missed;Entries", mBitPatternsAfter.size(), 0, mBitPatternsAfter.size(), mBitPatternsAfter.size(), 0, mBitPatternsAfter.size()); + makePatternAxisLabels(mHExtensionPatternsGoodMissed.get(), false); + mHExtensionPatternsGoodEmpty = createHistogram("hExtensionPatternsGoodEmpty", "Extended Tracks Pattern (GOOD) Empty Clusters;Before;After;Entries", mBitPatternsAfter.size(), 0, mBitPatternsAfter.size(), mBitPatternsAfter.size(), 0, mBitPatternsAfter.size()); + makePatternAxisLabels(mHExtensionPatternsGoodEmpty.get(), false); + + /// Effiencies + mEExtensionNum = createHistogram("hExtensionNum", "Extension Numerator", effHistBins, xbins.data()); + mEExtensionDen = createHistogram("hExtensionDen", "Extension Dennominator", effHistBins, xbins.data()); + // Purity + mEExtensionPurityNum = createHistogram("hExtensionPurityNum", "Extension Purity Numerator", effHistBins, xbins.data()); + mEExtensionPurityDen = createHistogram("hExtensionPurityDen", "Extension Purity Denominator", effHistBins, xbins.data()); + // Fake + mEExtensionFakeNum = createHistogram("hExtensionFakeNum", "Extension Fake Numerator", effHistBins, xbins.data()); + mEExtensionFakeDen = createHistogram("hExtensionFakeDen", "Extension Fake Denominator", effHistBins, xbins.data()); + mEExtensionFakeBeforeNum = createHistogram("hExtensionFakeBeforeNum", "Extension Fake Before Numerator", effHistBins, xbins.data()); + mEExtensionFakeAfterNum = createHistogram("hExtensionFakeAfterNum", "Extension Fake After Numerator", effHistBins, xbins.data()); + mEExtensionFakeMixNum = createHistogram("hExtensionFakeMixNum", "Extension Fake Mix Numerator", effHistBins, xbins.data()); + // Top + mEExtensionTopNum = createHistogram("hExtensionTopNum", "Extension Top Numerator", effHistBins, xbins.data()); + mEExtensionTopPurityNum = createHistogram("hExtensionTopPurityNum", "Extension Top Purity Numerator", effHistBins, xbins.data()); + mEExtensionTopFakeNum = createHistogram("hExtensionTopFakeNum", "Extension Top Fake Numerator", effHistBins, xbins.data()); + mEExtensionBotNum = createHistogram("hExtensionBotNum", "Extension Bot Numerator", effHistBins, xbins.data()); + mEExtensionBotPurityNum = createHistogram("hExtensionBotPurityNum", "Extension Bot Purity Numerator", effHistBins, xbins.data()); + mEExtensionBotFakeNum = createHistogram("hExtensionBotFakeNum", "Extension Bot Fake Numerator", effHistBins, xbins.data()); + mEExtensionMixNum = createHistogram("hExtensionMixNum", "Extension Mix Numerator", effHistBins, xbins.data()); + mEExtensionMixPurityNum = createHistogram("hExtensionMixPurityNum", "Extension Mix Purity Numerator", effHistBins, xbins.data()); + mEExtensionMixFakeNum = createHistogram("hExtensionMixFakeNum", "Extension Mix Fake Numerator", effHistBins, xbins.data()); + // Patterns + for (int i{0}; i < mBitPatternsBefore.size(); ++i) { + mEExtensionPatternGoodNum[i] = createHistogram(fmt::format("hExtensionPatternGood_{:07b}", mBitPatternsBefore[i]).c_str(), fmt::format("Extended Tracks Pattern (GOOD) {:07b}", mBitPatternsBefore[i]).c_str(), effHistBins, xbins.data()); + mEExtensionPatternFakeNum[i] = createHistogram(fmt::format("hExtensionPatternFake_{:07b}", mBitPatternsBefore[i]).c_str(), fmt::format("Extended Tracks Pattern (FAKE) {:07b}", mBitPatternsBefore[i]).c_str(), effHistBins, xbins.data()); + for (int j{0}; j < mBitPatternsAfter.size(); ++j) { + mEExtensionPatternIndGoodNum[i][j] = createHistogram(fmt::format("hExtensionPatternGood_{:07b}_{:07b}", mBitPatternsBefore[i], mBitPatternsAfter[j]).c_str(), fmt::format("Extended Tracks Pattern (GOOD) {:07b} -> {:07b}", mBitPatternsBefore[i], mBitPatternsAfter[j]).c_str(), effHistBins, xbins.data()); + mEExtensionPatternIndFakeNum[i][j] = createHistogram(fmt::format("hExtensionPatternFake_{:07b}_{:07b}", mBitPatternsBefore[i], mBitPatternsAfter[j]).c_str(), fmt::format("Extended Tracks Pattern (FAKE) {:07b} -> {:07b}", mBitPatternsBefore[i], mBitPatternsAfter[j]).c_str(), effHistBins, xbins.data()); + } + } + + /// DCA + mDCAxyVsPtPionsNormal = createHistogram("hDCAxyVsPtResNormal", "DCA_{#it{xy}} NORMAL Pions;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", effHistBins, xbins.data(), 1000, -500, 500); + mDCAxyVsPtPionsExtended = createHistogram("hDCAxyVsPtResExtended", "DCA_{#it{xy}} EXTENDED Pions;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", effHistBins, xbins.data(), 1000, -500, 500); + mDCAzVsPtPionsNormal = createHistogram("hDCAzVsPtResNormal", "DCA_{#it{z}} NORMAL Pions;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", effHistBins, xbins.data(), 1000, -500, 500); + mDCAzVsPtPionsExtended = createHistogram("hDCAzVsPtResExtended", "DCA_{#it{z}} EXTENDED Pions;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", effHistBins, xbins.data(), 1000, -500, 500); mStream = std::make_unique(mOutFileName.c_str(), "RECREATE"); } @@ -192,8 +257,6 @@ void TrackExtensionStudy::run(ProcessingContext& pc) mClustersMCLCont = recoData.getITSClustersMCLabels(); mInputITSidxs = recoData.getITSTracksClusterRefs(); - mHClustersCounts->SetBinContent(1, mClusters.size()); - LOGP(info, "** Found in {} rofs:\n\t- {} clusters with {} labels\n\t- {} tracks with {} labels", mTracksROFRecords.size(), mClusters.size(), mClustersMCLCont->getIndexedSize(), mTracks.size(), mTracksMCLabels.size()); LOGP(info, "** Found {} sources from kinematic files", mKineReader->getNSources()); @@ -208,10 +271,13 @@ void TrackExtensionStudy::process() for (int iSource{0}; iSource < mKineReader->getNSources(); ++iSource) { mParticleInfo[iSource].resize(mKineReader->getNEvents(iSource)); // events for (int iEvent{0}; iEvent < mKineReader->getNEvents(iSource); ++iEvent) { + const auto& mcEvent = mKineReader->getMCEventHeader(iSource, iEvent); mParticleInfo[iSource][iEvent].resize(mKineReader->getTracks(iSource, iEvent).size()); // tracks for (auto iPart{0}; iPart < mKineReader->getTracks(iEvent).size(); ++iPart) { - auto& part = mKineReader->getTracks(iSource, iEvent)[iPart]; - mParticleInfo[iSource][iEvent][iPart].event = iEvent; + const auto& part = mKineReader->getTracks(iSource, iEvent)[iPart]; + mParticleInfo[iSource][iEvent][iPart].eventX = mcEvent.GetX(); + mParticleInfo[iSource][iEvent][iPart].eventY = mcEvent.GetY(); + mParticleInfo[iSource][iEvent][iPart].eventZ = mcEvent.GetZ(); mParticleInfo[iSource][iEvent][iPart].pdg = part.GetPdgCode(); mParticleInfo[iSource][iEvent][iPart].pt = part.GetPt(); mParticleInfo[iSource][iEvent][iPart].phi = part.GetPhi(); @@ -287,6 +353,8 @@ void TrackExtensionStudy::process() LOGP(info, "\t- Total number of good: {} ({:.2f} %)", good, good * 100. / mTracks.size()); LOGP(info, "\t- Total number of extensions: {} ({:.2f} %)", extended, extended * 100. / mTracks.size()); + o2::dataformats::VertexBase collision; + o2::dataformats::DCA impactParameter; LOGP(info, "** Filling histograms ... "); for (auto iTrack{0}; iTrack < mTracks.size(); ++iTrack) { auto& lab = mTracksMCLabels[iTrack]; @@ -302,6 +370,7 @@ void TrackExtensionStudy::process() continue; } const auto& trk = part.track; + bool isGood = part.isReco && !part.isFake; mHTrackCounts->Fill(0); std::bitset<7> extPattern{0}; @@ -311,44 +380,85 @@ void TrackExtensionStudy::process() } } - if (!extPattern.any()) { - mHTrackCounts->Fill(1); - if (part.isReco || !part.isFake) { - mHTrackCounts->Fill(2); - } else { - mHTrackCounts->Fill(3); + // Tree + while (mWithTree) { + constexpr float refRadius{70.f}; + constexpr float maxSnp{0.9f}; + auto cTrk = trk; + if (!o2::base::Propagator::Instance()->PropagateToXBxByBz(cTrk, refRadius, maxSnp, 2.f, o2::base::Propagator::MatCorrType::USEMatCorrTGeo)) { + break; } - continue; - } - - if (mWithTree) { std::array xyz{(float)part.mcTrack.GetStartVertexCoordinatesX(), (float)part.mcTrack.GetStartVertexCoordinatesY(), (float)part.mcTrack.GetStartVertexCoordinatesZ()}; std::array pxyz{(float)part.mcTrack.GetStartVertexMomentumX(), (float)part.mcTrack.GetStartVertexMomentumY(), (float)part.mcTrack.GetStartVertexMomentumZ()}; auto pdg = O2DatabasePDG::Instance()->GetParticle(part.pdg); if (pdg == nullptr) { - LOGP(fatal, "MC info not available"); + LOGP(error, "MC info not available"); + break; } auto mcTrk = o2::track::TrackPar(xyz, pxyz, TMath::Nint(pdg->Charge() / 3.), true); + if (!mcTrk.rotate(cTrk.getAlpha()) || !o2::base::Propagator::Instance()->PropagateToXBxByBz(mcTrk, refRadius, maxSnp, 2.f, o2::base::Propagator::MatCorrType::USEMatCorrTGeo)) { + break; + } (*mStream) << "tree" - << "trk=" << trk + << "trk=" << cTrk << "mcTrk=" << mcTrk + << "isGood=" << isGood + << "isExtended=" << extPattern.any() << "\n"; + break; + } + + // impact parameter + while (isGood && std::abs(part.pdg) == 211) { + auto trkC = part.track; + collision.setXYZ(part.eventX, part.eventY, part.eventZ); + if (!o2::base::Propagator::Instance()->propagateToDCA(collision, trkC, o2::base::Propagator::Instance()->getNominalBz(), 2.0, o2::base::Propagator::MatCorrType::USEMatCorrTGeo, &impactParameter)) { + break; + } + + auto dcaXY = impactParameter.getY() * 1e4; + auto dcaZ = impactParameter.getZ() * 1e4; + if (!extPattern.any()) { + mDCAxyVsPtPionsNormal->Fill(part.pt, dcaXY); + mDCAzVsPtPionsNormal->Fill(part.pt, dcaZ); + } else { + mDCAxyVsPtPionsExtended->Fill(part.pt, dcaXY); + mDCAzVsPtPionsExtended->Fill(part.pt, dcaZ); + } + break; + } + + mEExtensionDen->Fill(trk.getPt()); + + if (!extPattern.any()) { + mHTrackCounts->Fill(1); + if (part.isReco || !part.isFake) { + mHTrackCounts->Fill(2); + } else { + mHTrackCounts->Fill(3); + } + continue; } mHTrackCounts->Fill(4); mHLengthAny->Fill(trk.getNClusters()); mHChi2Any->Fill(trk.getChi2()); mHPtAny->Fill(trk.getPt()); - if (part.isReco || !part.isFake) { + mEExtensionNum->Fill(trk.getPt()); + mEExtensionPurityDen->Fill(trk.getPt()); + mEExtensionFakeDen->Fill(trk.getPt()); + if (isGood) { mHTrackCounts->Fill(5); mHLengthGood->Fill(trk.getNClusters()); mHChi2Good->Fill(trk.getChi2()); mHPtGood->Fill(trk.getPt()); + mEExtensionPurityNum->Fill(trk.getPt()); } else { mHTrackCounts->Fill(6); mHLengthFake->Fill(trk.getNClusters()); mHChi2Fake->Fill(trk.getChi2()); mHPtFake->Fill(trk.getPt()); + mEExtensionFakeNum->Fill(trk.getPt()); } std::bitset<7> clusPattern{static_cast(trk.getPattern())}; @@ -356,29 +466,28 @@ void TrackExtensionStudy::process() if (extPattern.test(iLayer)) { extPattern.set(iLayer); mHExtensionAny->Fill(iLayer); - if (part.isReco || !part.isFake) { + if (isGood) { mHExtensionGood->Fill(iLayer); } else { mHExtensionFake->Fill(iLayer); } - - if ((part.fakeClusters & (0x1 << iLayer)) == 0) { - mHClustersCounts->Fill(4); - } else { - mHClustersCounts->Fill(5); - } } } - std::bitset<7> oldPattern{clusPattern & ~extPattern}; + std::bitset<7> oldPattern{clusPattern & ~extPattern}, holePattern{clusPattern}; + holePattern.flip(); auto clusN = clusPattern.to_ulong(); auto clusIdx = std::distance(std::begin(mBitPatternsAfter), std::find(std::begin(mBitPatternsAfter), std::end(mBitPatternsAfter), clusN)); auto oldN = oldPattern.to_ulong(); auto oldIdx = std::distance(std::begin(mBitPatternsBefore), std::find(std::begin(mBitPatternsBefore), std::end(mBitPatternsBefore), oldN)); mHExtensionPatternsAny->Fill(oldIdx, clusIdx); - if (part.isReco || !part.isFake) { + if (isGood) { mHExtensionPatternsGood->Fill(oldIdx, clusIdx); + mEExtensionPatternGoodNum[oldIdx]->Fill(trk.getPt()); + mEExtensionPatternIndGoodNum[oldIdx][clusIdx]->Fill(trk.getPt()); } else { mHExtensionPatternsFake->Fill(oldIdx, clusIdx); + mEExtensionPatternFakeNum[oldIdx]->Fill(trk.getPt()); + mEExtensionPatternIndFakeNum[oldIdx][clusIdx]->Fill(trk.getPt()); } // old pattern @@ -392,17 +501,70 @@ void TrackExtensionStudy::process() } } } - if (oldFake && newFake) { mHTrackCounts->Fill(9); + mEExtensionFakeMixNum->Fill(trk.getPt()); } else if (oldFake) { mHTrackCounts->Fill(7); + mEExtensionFakeBeforeNum->Fill(trk.getPt()); } else if (newFake) { mHTrackCounts->Fill(8); + mEExtensionFakeAfterNum->Fill(trk.getPt()); } - mHClustersCounts->SetBinContent(2, mHClustersCounts->GetBinContent(2) + oldPattern.count()); - mHClustersCounts->SetBinContent(3, mHClustersCounts->GetBinContent(3) + extPattern.count()); + // Check if we missed some clusters + if (isGood && holePattern.any()) { + auto missPattern{clusPattern}, emptyPattern{clusPattern}; + for (int iLayer{0}; iLayer < 7; ++iLayer) { + if (!holePattern.test(iLayer)) { + continue; + } + + // Check if there was actually a cluster that we missed + if ((part.clusters & (1 << iLayer)) != 0) { + missPattern.set(iLayer); + } else { + emptyPattern.set(iLayer); + } + } + + if (missPattern != clusPattern) { + auto missN = missPattern.to_ulong(); + auto missIdx = std::distance(std::begin(mBitPatternsAfter), std::find(std::begin(mBitPatternsAfter), std::end(mBitPatternsAfter), missN)); + mHExtensionPatternsGoodMissed->Fill(clusIdx, missIdx); + } + if (emptyPattern != clusPattern) { + auto emptyN = emptyPattern.to_ulong(); + auto emptyIdx = std::distance(std::begin(mBitPatternsAfter), std::find(std::begin(mBitPatternsAfter), std::end(mBitPatternsAfter), emptyN)); + mHExtensionPatternsGoodEmpty->Fill(clusIdx, emptyIdx); + } + } + + // Top/Bot/Mixed Extension + bool isTop = (extPattern & mTopMask).any(); + bool isBot = (extPattern & mBotMask).any(); + if (isTop && isBot) { + mEExtensionMixNum->Fill(trk.getPt()); + if (isGood) { + mEExtensionMixPurityNum->Fill(trk.getPt()); + } else { + mEExtensionMixFakeNum->Fill(trk.getPt()); + } + } else if (isTop) { + mEExtensionTopNum->Fill(trk.getPt()); + if (isGood) { + mEExtensionTopPurityNum->Fill(trk.getPt()); + } else { + mEExtensionTopFakeNum->Fill(trk.getPt()); + } + } else { + mEExtensionBotNum->Fill(trk.getPt()); + if (isGood) { + mEExtensionBotPurityNum->Fill(trk.getPt()); + } else { + mEExtensionBotFakeNum->Fill(trk.getPt()); + } + } } } @@ -421,39 +583,57 @@ void TrackExtensionStudy::endOfStream(EndOfStreamContext& ec) { LOGP(info, "Writing results to {}", mOutFileName); mStream->GetFile()->cd(); + for (const auto h : mHistograms) { + h->Write(); + } - mHTrackCounts->Write(); - mHClustersCounts->Write(); - - mHLengthAny->Write(); - mHLengthGood->Write(); - mHLengthFake->Write(); - - mHChi2Any->Write(); - mHChi2Good->Write(); - mHChi2Fake->Write(); - - mHPtAny->Write(); - mHPtGood->Write(); - mHPtFake->Write(); - - mHExtensionAny->Write(); - mHExtensionGood->Write(); - mHExtensionFake->Write(); - - mHExtensionPatternsAny->Write(); - mHExtensionPatternsGood->Write(); - mHExtensionPatternsFake->Write(); + LOGP(info, "Calculating efficiencies"); + auto makeEff = [](auto num, auto den, const char* name, const char* title) { + auto e = std::make_unique(*num, *den); + e->SetName(name); + e->SetTitle(title); + e->Write(); + }; + makeEff(mEExtensionNum.get(), mEExtensionDen.get(), "eExtension", "Track Extension EXT TRK/ALL"); + makeEff(mEExtensionPurityNum.get(), mEExtensionPurityDen.get(), "eExtensionPurity", "Track Extension Purity GOOD/EXT TRK"); + makeEff(mEExtensionFakeNum.get(), mEExtensionFakeDen.get(), "eExtensionFake", "Track Extension Fake FAKE/EXT TRK"); + makeEff(mEExtensionFakeBeforeNum.get(), mEExtensionFakeNum.get(), "eExtensionFakeBefore", "Track Extension Fake FAKE BEF/FAKE EXT TRK"); + makeEff(mEExtensionFakeAfterNum.get(), mEExtensionFakeNum.get(), "eExtensionFakeAfter", "Track Extension Fake FAKE AFT/FAKE EXT TRK"); + makeEff(mEExtensionFakeMixNum.get(), mEExtensionFakeNum.get(), "eExtensionFakeMix", "Track Extension Fake FAKE MIX/FAKE EXT TRK"); + makeEff(mEExtensionTopNum.get(), mEExtensionDen.get(), "eExtensionTop", "Track Extension Top"); + makeEff(mEExtensionTopPurityNum.get(), mEExtensionPurityDen.get(), "eExtensionTopPurity", "Track Extension Purity GOOD TOP/EXT TRK"); + makeEff(mEExtensionTopFakeNum.get(), mEExtensionFakeNum.get(), "eExtensionTopFake", "Track Extension FAKE TOP/EXT FAKE TRK"); + makeEff(mEExtensionBotNum.get(), mEExtensionDen.get(), "eExtensionBot", "Track Extension Bot"); + makeEff(mEExtensionBotPurityNum.get(), mEExtensionPurityDen.get(), "eExtensionBotPurity", "Track Extension Purity GOOD BOT/EXT TRK"); + makeEff(mEExtensionBotFakeNum.get(), mEExtensionFakeNum.get(), "eExtensionBotFake", "Track Extension FAKE BOT/EXT FAKE TRK"); + makeEff(mEExtensionMixNum.get(), mEExtensionDen.get(), "eExtensionMix", "Track Extension Mix"); + makeEff(mEExtensionMixPurityNum.get(), mEExtensionPurityDen.get(), "eExtensionMixPurity", "Track Extension Purity GOOD MIX/EXT TRK"); + makeEff(mEExtensionMixFakeNum.get(), mEExtensionFakeNum.get(), "eExtensionMixFake", "Track Extension FAKE MIX/EXT FAKE TRK"); + for (int i{0}; i < mBitPatternsBefore.size(); ++i) { + makeEff(mEExtensionPatternGoodNum[i].get(), mEExtensionPurityNum.get(), fmt::format("eExtensionPatternGood_{:07b}", mBitPatternsBefore[i]).c_str(), fmt::format("Extended Tracks Pattern (GOOD) {:07b} GOOD EXT TRK/EXT TRK", mBitPatternsBefore[i]).c_str()); + makeEff(mEExtensionPatternFakeNum[i].get(), mEExtensionFakeNum.get(), fmt::format("eExtensionPatternFake_{:07b}", mBitPatternsBefore[i]).c_str(), fmt::format("Extended Tracks Pattern (FAKE) {:07b} FAKE EXT TRK/EXT TRK", mBitPatternsBefore[i]).c_str()); + for (int j{0}; j < mBitPatternsAfter.size(); ++j) { + makeEff(mEExtensionPatternIndGoodNum[i][j].get(), mEExtensionPatternGoodNum[i].get(), fmt::format("eExtensionPatternGood_{:07b}_{:07b}", mBitPatternsBefore[i], mBitPatternsAfter[j]).c_str(), fmt::format("Extended Tracks Pattern (GOOD) {:07b} -> {:07b} GOOD EXT TRK/EXT TRK", mBitPatternsBefore[i], mBitPatternsAfter[j]).c_str()); + makeEff(mEExtensionPatternIndFakeNum[i][j].get(), mEExtensionPatternFakeNum[i].get(), fmt::format("eExtensionPatternFake_{:07b}_{:07b}", mBitPatternsBefore[i], mBitPatternsAfter[j]).c_str(), fmt::format("Extended Tracks Pattern (FAKE) {:07b} -> {:07b} FAKE EXT TRK/EXT TRK", mBitPatternsBefore[i], mBitPatternsAfter[j]).c_str()); + } + } mStream->Close(); } -DataProcessorSpec getTrackExtensionStudy(mask_t srcTracksMask, mask_t srcClustersMask, bool useMC, std::shared_ptr kineReader) +void TrackExtensionStudy::finaliseCCDB(ConcreteDataMatcher& matcher, void* obj) +{ + if (o2::base::GRPGeomHelper::instance().finaliseCCDB(matcher, obj)) { + return; + } +} + +DataProcessorSpec getTrackExtensionStudy(mask_t srcTracksMask, mask_t srcClustersMask, std::shared_ptr kineReader) { std::vector outputs; auto dataRequest = std::make_shared(); - dataRequest->requestTracks(srcTracksMask, useMC); - dataRequest->requestClusters(srcClustersMask, useMC); + dataRequest->requestTracks(srcTracksMask, true); + dataRequest->requestClusters(srcClustersMask, true); auto ggRequest = std::make_shared(false, // orbitResetTime true, // GRPECS=true @@ -468,7 +648,7 @@ DataProcessorSpec getTrackExtensionStudy(mask_t srcTracksMask, mask_t srcCluster "its-study-track-extension", dataRequest->inputs, outputs, - AlgorithmSpec{adaptFromTask(dataRequest, srcTracksMask, useMC, kineReader, ggRequest)}, + AlgorithmSpec{adaptFromTask(dataRequest, srcTracksMask, kineReader, ggRequest)}, Options{{"with-tree", o2::framework::VariantType::Bool, false, {"Produce in addition a tree"}}}}; } diff --git a/Detectors/ITSMFT/ITS/postprocessing/workflow/standalone-postprocessing-workflow.cxx b/Detectors/ITSMFT/ITS/postprocessing/workflow/standalone-postprocessing-workflow.cxx index 02a75def154fc..30fb39c77f235 100644 --- a/Detectors/ITSMFT/ITS/postprocessing/workflow/standalone-postprocessing-workflow.cxx +++ b/Detectors/ITSMFT/ITS/postprocessing/workflow/standalone-postprocessing-workflow.cxx @@ -113,11 +113,14 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) specs.emplace_back(o2::its::study::getAnomalyStudy(srcCls, useMC)); } if (configcontext.options().get("track-extension-study")) { + if (!useMC) { + LOGP(fatal, "Track Extension Study needs MC!"); + } anyStudy = true; srcTrc = GID::getSourcesMask(configcontext.options().get("track-sources")); srcCls = GID::getSourcesMask("ITS"); - o2::globaltracking::InputHelper::addInputSpecs(configcontext, specs, srcCls, srcTrc, srcTrc, useMC, srcCls, srcTrc); - specs.emplace_back(o2::its::study::getTrackExtensionStudy(srcTrc, srcCls, useMC, mcKinematicsReader)); + o2::globaltracking::InputHelper::addInputSpecs(configcontext, specs, srcCls, srcTrc, srcTrc, true, srcCls, srcTrc); + specs.emplace_back(o2::its::study::getTrackExtensionStudy(srcTrc, srcCls, mcKinematicsReader)); } if (configcontext.options().get("efficiency-study")) { anyStudy = true; diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h index 1c4d3360bc7a3..976d01f1d476b 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h @@ -94,11 +94,17 @@ struct TrackingParameters { unsigned long MaxMemory = 12000000000UL; float MaxChi2ClusterAttachment = 60.f; float MaxChi2NDF = 30.f; - bool UseTrackFollower = false; bool FindShortTracks = false; bool PerPrimaryVertexProcessing = false; bool SaveTimeBenchmarks = false; bool DoUPCIteration = false; + /// Cluster attachment + bool UseTrackFollower = false; + bool UseTrackFollowerTop = false; + bool UseTrackFollowerBot = false; + bool UseTrackFollowerMix = false; + float TrackFollowerNSigmaCutZ = 1.f; + float TrackFollowerNSigmaCutPhi = 1.f; }; inline int TrackingParameters::CellMinimumLevel() diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracker.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracker.h index 70de43d83d8d2..58483e4aa9f6f 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracker.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracker.h @@ -122,7 +122,7 @@ float Tracker::evaluateTask(void (Tracker::*task)(T...), const char* taskName, s { float diff{0.f}; - if (constants::DoTimeBenchmarks) { + if constexpr (constants::DoTimeBenchmarks) { auto start = std::chrono::high_resolution_clock::now(); (this->*task)(std::forward(args)...); auto end = std::chrono::high_resolution_clock::now(); diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraits.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraits.h index 207dd5d3d50f5..46499db92d4d5 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraits.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraits.h @@ -78,10 +78,10 @@ class TrackerTraits bool isMatLUT() const; // Others - GPUhd() static constexpr int4 getEmptyBinsRect() { return int4{0, 0, 0, 0}; } - const int4 getBinsRect(const Cluster&, int layer, float z1, float z2, float maxdeltaz, float maxdeltaphi); - const int4 getBinsRect(int layer, float phi, float maxdeltaphi, float z, float maxdeltaz); - const int4 getBinsRect(int layer, float phi, float maxdeltaphi, float z1, float z2, float maxdeltaz); + GPUhd() static consteval int4 getEmptyBinsRect() { return int4{0, 0, 0, 0}; } + const int4 getBinsRect(const Cluster&, int layer, float z1, float z2, float maxdeltaz, float maxdeltaphi) const noexcept; + const int4 getBinsRect(int layer, float phi, float maxdeltaphi, float z, float maxdeltaz) const noexcept; + const int4 getBinsRect(int layer, float phi, float maxdeltaphi, float z1, float z2, float maxdeltaz) const noexcept; void SetRecoChain(o2::gpu::GPUChainITS* chain) { mChain = chain; } void setSmoothing(bool v) { mApplySmoothing = v; } bool getSmoothing() const { return mApplySmoothing; } @@ -112,6 +112,12 @@ class TrackerTraits bool mIsGPU = false; }; +inline void TrackerTraits::initialiseTimeFrame(const int iteration) +{ + mTimeFrame->initialise(iteration, mTrkParams[iteration], mTrkParams[iteration].NLayers); + setIsGPU(false); +} + inline float TrackerTraits::getBz() const { return mBz; @@ -122,40 +128,32 @@ inline void TrackerTraits::UpdateTrackingParameters(const std::vectorinitialise(iteration, mTrkParams[iteration], mTrkParams[iteration].NLayers); - setIsGPU(false); -} - -inline const int4 TrackerTraits::getBinsRect(const int layerIndex, float phi, float maxdeltaphi, - float z1, float z2, float maxdeltaz) +inline const int4 TrackerTraits::getBinsRect(const int layerIndex, float phi, float maxdeltaphi, float z1, float z2, float maxdeltaz) const noexcept { const float zRangeMin = o2::gpu::GPUCommonMath::Min(z1, z2) - maxdeltaz; const float phiRangeMin = (maxdeltaphi > constants::math::Pi) ? 0.f : phi - maxdeltaphi; const float zRangeMax = o2::gpu::GPUCommonMath::Max(z1, z2) + maxdeltaz; const float phiRangeMax = (maxdeltaphi > constants::math::Pi) ? constants::math::TwoPi : phi + maxdeltaphi; - if (zRangeMax < -mTrkParams[0].LayerZ[layerIndex + 1] || - zRangeMin > mTrkParams[0].LayerZ[layerIndex + 1] || zRangeMin > zRangeMax) { - + if (zRangeMax < -mTrkParams[0].LayerZ[layerIndex] || + zRangeMin > mTrkParams[0].LayerZ[layerIndex] || zRangeMin > zRangeMax) { return getEmptyBinsRect(); } const IndexTableUtils& utils{mTimeFrame->mIndexTableUtils}; - return int4{o2::gpu::GPUCommonMath::Max(0, utils.getZBinIndex(layerIndex + 1, zRangeMin)), + return int4{o2::gpu::GPUCommonMath::Max(0, utils.getZBinIndex(layerIndex, zRangeMin)), utils.getPhiBinIndex(math_utils::getNormalizedPhi(phiRangeMin)), - o2::gpu::GPUCommonMath::Min(mTrkParams[0].ZBins - 1, utils.getZBinIndex(layerIndex + 1, zRangeMax)), // /!\ trkParams can potentially change across iterations + o2::gpu::GPUCommonMath::Min(mTrkParams[0].ZBins - 1, utils.getZBinIndex(layerIndex, zRangeMax)), // /!\ trkParams can potentially change across iterations utils.getPhiBinIndex(math_utils::getNormalizedPhi(phiRangeMax))}; } } // namespace its diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingConfigParam.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingConfigParam.h index fe5e52bd6277a..68bfdb51170b5 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingConfigParam.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingConfigParam.h @@ -72,7 +72,9 @@ struct TrackerParamConfig : public o2::conf::ConfigurableParamHelper 0 off + float trackFollowerNSigmaZ = 1.f; // sigma in z-cut for track-following search rectangle + float trackFollowerNSigmaPhi = 1.f; // sigma in phi-cut for track-following search rectangle float cellsPerClusterLimit = -1.f; float trackletsPerClusterLimit = -1.f; int findShortTracks = -1; diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Vertexer.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Vertexer.h index 7fb5d1421877e..ac0cf51921176 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Vertexer.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Vertexer.h @@ -176,7 +176,7 @@ float Vertexer::evaluateTask(void (Vertexer::*task)(T...), const char* taskName, { float diff{0.f}; - if (constants::DoTimeBenchmarks) { + if constexpr (constants::DoTimeBenchmarks) { auto start = std::chrono::high_resolution_clock::now(); (this->*task)(std::forward(args)...); auto end = std::chrono::high_resolution_clock::now(); diff --git a/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx b/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx index 9a43402a2e93a..721452bf0361d 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx @@ -120,7 +120,7 @@ void Tracker::clustersToTracks(std::function logger, std::f total += evaluateTask(&Tracker::findShortPrimaries, "Short primaries finding", logger); std::stringstream sstream; - if (constants::DoTimeBenchmarks) { + if constexpr (constants::DoTimeBenchmarks) { sstream << std::setw(2) << " - " << "Timeframe " << mTimeFrameCounter++ << " processing completed in: " << total << "ms using " << mTraits->getNThreads() << " threads."; } @@ -200,7 +200,7 @@ void Tracker::clustersToTracksHybrid(std::function logger, // total += evaluateTask(&Tracker::findShortPrimaries, "Hybrid short primaries finding", logger); std::stringstream sstream; - if (constants::DoTimeBenchmarks) { + if constexpr (constants::DoTimeBenchmarks) { sstream << std::setw(2) << " - " << "Timeframe " << mTimeFrameCounter++ << " processing completed in: " << total << "ms using " << mTraits->getNThreads() << " threads."; } @@ -502,8 +502,16 @@ void Tracker::getGlobalConfiguration() if (tc.maxMemory) { params.MaxMemory = tc.maxMemory; } - if (tc.useTrackFollower >= 0) { - params.UseTrackFollower = tc.useTrackFollower; + if (tc.useTrackFollower > 0) { + params.UseTrackFollower = true; + // Bit 0: Allow for mixing of top&bot extension --> implies Bits 1&2 set + // Bit 1: Allow for top extension + // Bit 2: Allow for bot extension + params.UseTrackFollowerMix = ((tc.useTrackFollower & (1 << 0)) != 0); + params.UseTrackFollowerTop = ((tc.useTrackFollower & (1 << 1)) != 0); + params.UseTrackFollowerBot = ((tc.useTrackFollower & (1 << 2)) != 0); + params.TrackFollowerNSigmaCutZ = tc.trackFollowerNSigmaZ; + params.TrackFollowerNSigmaCutPhi = tc.trackFollowerNSigmaPhi; } if (tc.cellsPerClusterLimit >= 0) { params.CellsPerClusterLimit = tc.cellsPerClusterLimit; diff --git a/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx b/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx index 4457d4515e0a6..da0abbae9dc1f 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx @@ -111,7 +111,7 @@ void TrackerTraits::computeLayerTracklets(const int iteration, int iROFslice, in const float sqInverseDeltaZ0{1.f / (Sq(currentCluster.zCoordinate - primaryVertex.getZ()) + 2.e-8f)}; /// protecting from overflows adding the detector resolution const float sigmaZ{o2::gpu::CAMath::Sqrt(Sq(resolution) * Sq(tanLambda) * ((Sq(inverseR0) + sqInverseDeltaZ0) * Sq(meanDeltaR) + 1.f) + Sq(meanDeltaR * tf->getMSangle(iLayer)))}; - const int4 selectedBinsRect{getBinsRect(currentCluster, iLayer, zAtRmin, zAtRmax, + const int4 selectedBinsRect{getBinsRect(currentCluster, iLayer + 1, zAtRmin, zAtRmax, sigmaZ * mTrkParams[iteration].NSigmaCut, tf->getPhiCut(iLayer))}; if (selectedBinsRect.x == 0 && selectedBinsRect.y == 0 && selectedBinsRect.z == 0 && selectedBinsRect.w == 0) { continue; @@ -679,10 +679,11 @@ void TrackerTraits::extendTracks(const int iteration) for (auto& track : mTimeFrame->getTracks(rof)) { auto backup{track}; bool success{false}; - if (track.getLastClusterLayer() != mTrkParams[iteration].NLayers - 1) { + // the order here biases towards top extension, tracks should probably be fitted separately in the directions and then compared. + if ((mTrkParams[iteration].UseTrackFollowerMix || mTrkParams[iteration].UseTrackFollowerTop) && track.getLastClusterLayer() != mTrkParams[iteration].NLayers - 1) { success = success || trackFollowing(&track, rof, true, iteration); } - if (track.getFirstClusterLayer() != 0) { + if ((mTrkParams[iteration].UseTrackFollowerMix || (mTrkParams[iteration].UseTrackFollowerBot && !success)) && track.getFirstClusterLayer() != 0) { success = success || trackFollowing(&track, rof, false, iteration); } if (success) { @@ -830,8 +831,8 @@ bool TrackerTraits::fitTrack(TrackITSExt& track, int start, int end, int step, f } if (mCorrType == o2::base::PropagatorF::MatCorrType::USEMatCorrNONE) { - float radl = 9.36f; // Radiation length of Si [cm] - float rho = 2.33f; // Density of Si [g/cm^3] + constexpr float radl = 9.36f; // Radiation length of Si [cm] + constexpr float rho = 2.33f; // Density of Si [g/cm^3] if (!track.correctForMaterial(mTrkParams[0].LayerxX0[iLayer], mTrkParams[0].LayerxX0[iLayer] * radl * rho, true)) { continue; } @@ -855,36 +856,40 @@ bool TrackerTraits::trackFollowing(TrackITSExt* track, int rof, bool outward, co auto propInstance = o2::base::Propagator::Instance(); const int step = -1 + outward * 2; const int end = outward ? mTrkParams[iteration].NLayers - 1 : 0; - std::vector hypotheses(1, *track); - for (auto& hypo : hypotheses) { - int iLayer = outward ? track->getLastClusterLayer() : track->getFirstClusterLayer(); + std::vector hypotheses(1, *track); // possibly avoid reallocation + for (size_t iHypo{0}; iHypo < hypotheses.size(); ++iHypo) { + auto hypo{hypotheses[iHypo]}; + int iLayer = static_cast(outward ? hypo.getLastClusterLayer() : hypo.getFirstClusterLayer()); + // per layer we add new hypotheses while (iLayer != end) { - iLayer += step; + iLayer += step; // step through all layers until we reach the end, this allows for skipping on empty layers const float r = mTrkParams[iteration].LayerRadii[iLayer]; + // get an estimate of the trackinf-frame x for the next step float x{-999}; if (!hypo.getXatLabR(r, x, mTimeFrame->getBz(), o2::track::DirAuto) || x <= 0.f) { continue; } - + // estimate hypo's trk parameters at that x auto& hypoParam{outward ? hypo.getParamOut() : hypo.getParamIn()}; if (!propInstance->propagateToX(hypoParam, x, mTimeFrame->getBz(), PropagatorF::MAX_SIN_PHI, PropagatorF::MAX_STEP, mTrkParams[iteration].CorrType)) { continue; } - if (mTrkParams[iteration].CorrType == PropagatorF::MatCorrType::USEMatCorrNONE) { - float radl = 9.36f; // Radiation length of Si [cm] - float rho = 2.33f; // Density of Si [g/cm^3] + if (mTrkParams[iteration].CorrType == PropagatorF::MatCorrType::USEMatCorrNONE) { // account for material affects if propagator does not + constexpr float radl = 9.36f; // Radiation length of Si [cm] + constexpr float rho = 2.33f; // Density of Si [g/cm^3] if (!hypoParam.correctForMaterial(mTrkParams[iteration].LayerxX0[iLayer], mTrkParams[iteration].LayerxX0[iLayer] * radl * rho, true)) { continue; } } + + // calculate the search window on this layer const float phi{hypoParam.getPhi()}; const float ePhi{o2::gpu::CAMath::Sqrt(hypoParam.getSigmaSnp2() / hypoParam.getCsp2())}; const float z{hypoParam.getZ()}; const float eZ{o2::gpu::CAMath::Sqrt(hypoParam.getSigmaZ2())}; - const int4 selectedBinsRect{getBinsRect(iLayer, phi, mTrkParams[iteration].NSigmaCut * ePhi, z, mTrkParams[iteration].NSigmaCut * eZ)}; - + const int4 selectedBinsRect{getBinsRect(iLayer, phi, mTrkParams[iteration].TrackFollowerNSigmaCutPhi * ePhi, z, mTrkParams[iteration].TrackFollowerNSigmaCutZ * eZ)}; if (selectedBinsRect.x == 0 && selectedBinsRect.y == 0 && selectedBinsRect.z == 0 && selectedBinsRect.w == 0) { continue; } @@ -900,9 +905,8 @@ bool TrackerTraits::trackFollowing(TrackITSExt* track, int rof, bool outward, co continue; } - TrackITSExt currentHypo{hypo}, newHypo{hypo}; - bool first{true}; - for (int iPhiCount{0}; iPhiCount < phiBinsNum; iPhiCount++) { + // check all clusters in search windows for possible new hypotheses + for (int iPhiCount = 0; iPhiCount < phiBinsNum; iPhiCount++) { int iPhiBin = (selectedBinsRect.y + iPhiCount) % mTrkParams[iteration].PhiBins; const int firstBinIndex{mTimeFrame->mIndexTableUtils.getBinIndex(selectedBinsRect.x, iPhiBin)}; const int maxBinIndex{firstBinIndex + selectedBinsRect.z - selectedBinsRect.x + 1}; @@ -921,7 +925,7 @@ bool TrackerTraits::trackFollowing(TrackITSExt* track, int rof, bool outward, co const TrackingFrameInfo& trackingHit = mTimeFrame->getTrackingFrameInfoOnLayer(iLayer).at(nextCluster.clusterId); - TrackITSExt& tbupdated = first ? hypo : newHypo; + auto tbupdated{hypo}; auto& tbuParams = outward ? tbupdated.getParamOut() : tbupdated.getParamIn(); if (!tbuParams.rotate(trackingHit.alphaTrackingFrame)) { continue; @@ -942,12 +946,7 @@ bool TrackerTraits::trackFollowing(TrackITSExt* track, int rof, bool outward, co } tbupdated.setChi2(tbupdated.getChi2() + predChi2); /// This is wrong for outward propagation as the chi2 refers to inward parameters tbupdated.setExternalClusterIndex(iLayer, nextCluster.clusterId, true); - - if (!first) { - hypotheses.emplace_back(tbupdated); - newHypo = currentHypo; - } - first = false; + hypotheses.emplace_back(tbupdated); } } } diff --git a/Detectors/MUON/MCH/Base/include/MCHBase/ResponseParam.h b/Detectors/MUON/MCH/Base/include/MCHBase/ResponseParam.h index a545bb670d59b..6dba967694026 100644 --- a/Detectors/MUON/MCH/Base/include/MCHBase/ResponseParam.h +++ b/Detectors/MUON/MCH/Base/include/MCHBase/ResponseParam.h @@ -30,11 +30,11 @@ struct ResponseParam : public o2::conf::ConfigurableParamHelper { float pitchSt1 = 0.21; ///< anode-cathode pitch (cm) for station 1 float pitchSt2345 = 0.25; ///< anode-cathode pitch (cm) for station 2 to 5 - float mathiesonSqrtKx3St1 = 0.7000; ///< Mathieson parameter sqrt(K3) in x direction for station 1 - float mathiesonSqrtKx3St2345 = 0.7131; ///< Mathieson parameter sqrt(K3) in x direction for station 2 to 5 + float mathiesonSqrtKx3St1 = 0.5477; ///< Mathieson parameter sqrt(K3) in x direction for station 1 + float mathiesonSqrtKx3St2345 = 0.5477; ///< Mathieson parameter sqrt(K3) in x direction for station 2 to 5 - float mathiesonSqrtKy3St1 = 0.7550; ///< Mathieson parameter sqrt(K3) in y direction for station 1 - float mathiesonSqrtKy3St2345 = 0.7642; ///< Mathieson parameter sqrt(K3) in y direction for station 2 to 5 + float mathiesonSqrtKy3St1 = 0.5477; ///< Mathieson parameter sqrt(K3) in y direction for station 1 + float mathiesonSqrtKy3St2345 = 0.5477; ///< Mathieson parameter sqrt(K3) in y direction for station 2 to 5 float chargeSlopeSt1 = 25.; ///< charge slope used in E to charge conversion for station 1 float chargeSlopeSt2345 = 10.; ///< charge slope used in E to charge conversion for station 2 to 5 diff --git a/Detectors/MUON/MCH/DevIO/Digits/digits-sampler-workflow.cxx b/Detectors/MUON/MCH/DevIO/Digits/digits-sampler-workflow.cxx index 7f3819f110ba3..0184e1c78c0c6 100644 --- a/Detectors/MUON/MCH/DevIO/Digits/digits-sampler-workflow.cxx +++ b/Detectors/MUON/MCH/DevIO/Digits/digits-sampler-workflow.cxx @@ -27,6 +27,7 @@ #include #include #include +#include using namespace o2::framework; @@ -63,61 +64,64 @@ class DigitSamplerTask : public io::DigitIOBaseTask void outputAndClear(DataAllocator& out) { - printSummary(mDigits, mROFs, "-> to output"); + LOGP(info, "Sending {} rofs with {} digits", mROFs.size(), mDigits.size()); out.snapshot(OutputRef{"rofs"}, mROFs); out.snapshot(OutputRef{"digits"}, mDigits); mDigits.clear(); mROFs.clear(); } - bool shouldEnd() const + bool shouldEnd() { bool maxTFreached = mNofProcessedTFs >= mMaxNofTimeFrames; bool maxROFreached = mNofProcessedROFs >= mMaxNofROFs; - return !mReadIsOk || maxTFreached || maxROFreached; + bool lastTF = mInput.peek() == EOF; + return !mReadIsOk || lastTF || maxTFreached || maxROFreached; } void run(ProcessingContext& pc) { if (shouldEnd()) { - // output remaining data if any - if (mROFs.size() > 0) { - --mTFid; - outputAndClear(pc.outputs()); - } - pc.services().get().endOfStream(); - return; + throw std::invalid_argument("process should have ended already"); } std::vector rofs; std::vector digits; - mReadIsOk = mDigitSampler->read(digits, rofs); - if (!mReadIsOk) { - return; - } + while ((mReadIsOk = mDigitSampler->read(digits, rofs))) { + + // process the current input TF if requested + if (shouldProcess()) { + incNofProcessedTFs(); + mNofProcessedROFs += rofs.size(); + // append rofs to mROFs, but shift the indices by the amount of digits + // we have read so far. + auto offset = mDigits.size(); + std::transform(rofs.begin(), rofs.end(), std::back_inserter(mROFs), + [offset](ROFRecord r) { + r.setDataRef(r.getFirstIdx() + offset, r.getNEntries()); + return r; + }); + mDigits.insert(mDigits.end(), digits.begin(), digits.end()); + printSummary(mDigits, mROFs); + printFull(mDigits, mROFs); + } - if (shouldProcess()) { - incNofProcessedTFs(); - mNofProcessedROFs += rofs.size(); - // append rofs to mROFs, but shift the indices by the amount of digits - // we have read so far. - auto offset = mDigits.size(); - std::transform(rofs.begin(), rofs.end(), std::back_inserter(mROFs), - [offset](ROFRecord r) { - r.setDataRef(r.getFirstIdx() + offset, r.getNEntries()); - return r; - }); - mDigits.insert(mDigits.end(), digits.begin(), digits.end()); - printSummary(mDigits, mROFs); - printFull(mDigits, mROFs); - } + // increment the input TF id for the next one + incTFid(); - // output if we've accumulated enough ROFs - if (mROFs.size() >= mMinNumberOfROFsPerTF) { - outputAndClear(pc.outputs()); + // stop here if we've accumulated enough ROFs or TFs + if (mROFs.size() >= mMinNumberOfROFsPerTF || shouldEnd()) { + break; + } } - incTFid(); + // output whatever has been accumulated, even if empty + outputAndClear(pc.outputs()); + + if (shouldEnd()) { + pc.services().get().endOfStream(); + pc.services().get().readyToQuit(QuitRequest::Me); + } } }; diff --git a/Detectors/MUON/MID/Calibration/macros/build_rejectlist.C b/Detectors/MUON/MID/Calibration/macros/build_rejectlist.C index 7a395d2c099da..48391b4460687 100644 --- a/Detectors/MUON/MID/Calibration/macros/build_rejectlist.C +++ b/Detectors/MUON/MID/Calibration/macros/build_rejectlist.C @@ -50,6 +50,17 @@ struct RejectListStruct { std::vector rejectList{}; /// Bad channels }; +/// @brief Useful metadata +struct MDStruct { + long start = 0; /// Start validity + long end = 0; /// End validity + int runNumber = 0; /// Run number + std::string runType; /// Run Type + + bool operator<(const MDStruct& other) const { return start < other.start; } + bool operator==(const MDStruct& other) const { return start == other.start; } +}; + /// @brief Get timestamp in milliseconds /// @param timestamp Input timestamp (in s or ms) /// @return Timestamp in ms @@ -96,23 +107,33 @@ std::string timeRangeToString(long start, long end) /// @param end Query objects created not after /// @param api CDB api /// @param path CDB path -/// @return Vector of start validity of each object sorted in ascending way -std::vector findObjectsTSInPeriod(long start, long end, const o2::ccdb::CcdbApi& api, const char* path) +/// @return Vector of metadata in ascending order +std::vector findObjectsMDInPeriod(long start, long end, const o2::ccdb::CcdbApi& api, const char* path) { - std::vector ts; - auto out = api.list(path, false, "text/plain", getTSMS(end), getTSMS(start)); - std::stringstream ss(out); - std::string token; - while (ss >> token) { - if (token.find("Validity") != std::string::npos) { - ss >> token; - ts.emplace_back(std::atol(token.c_str())); - } + std::vector mds; + auto out = api.list(path, false, "application/json", getTSMS(end), getTSMS(start)); + rapidjson::Document doc; + doc.Parse(out.c_str()); + for (auto& obj : doc["objects"].GetArray()) { + MDStruct md; + md.start = obj["validFrom"].GetInt64(); + md.end = obj["validUntil"].GetInt64(); + md.runNumber = std::atoi(obj["RunNumber"].GetString()); + md.runType = obj["RunType"].GetString(); + mds.emplace_back(md); } - ts.erase(std::unique(ts.begin(), ts.end()), ts.end()); + mds.erase(std::unique(mds.begin(), mds.end()), mds.end()); // Sort timestamps in ascending order - std::sort(ts.begin(), ts.end()); - return ts; + std::sort(mds.begin(), mds.end()); + return mds; +} + +/// @brief Gets the quality trend graph from the quality canvas +/// @param qcQuality MID QC quality canvas +/// @return Quality trend graph +TGraph* getQualityTrend(const TCanvas* qcQuality) +{ + return static_cast(qcQuality->GetListOfPrimitives()->FindObject("Graph")); } /// @brief Find the first and last time when the quality was good or bad @@ -127,7 +148,7 @@ std::pair findTSRange(TCanvas* qcQuality, bool selectBad = t // Medium: 2.5 // Bad: 1.5 // Null: 0.5 - auto* gr = static_cast(qcQuality->GetListOfPrimitives()->FindObject("Graph")); + auto* gr = getQualityTrend(qcQuality); double xp, yp; std::pair range{std::numeric_limits::max(), 0}; for (int ip = 0; ip < gr->GetN(); ++ip) { @@ -144,6 +165,32 @@ std::pair findTSRange(TCanvas* qcQuality, bool selectBad = t return range; } +/// @brief Gets the first and last timestamp in the quality +/// @param qcQuality MID QC quality canvas +/// @return Pair with the first and last timestamp in the quality trend +std::pair getFirstLast(const TCanvas* qcQuality) +{ + auto* gr = getQualityTrend(qcQuality); + double xp1, xp2, yp; + gr->GetPoint(0, xp1, yp); + gr->GetPoint(gr->GetN() - 1, xp2, yp); + return {static_cast(xp1 * 1000), static_cast(xp2 * 1000)}; +} + +/// @brief Update the selected range of timestamp +/// @param selectedTSRange Reference to the selected range to be modified +/// @param qcTSRange Range of the MID quality trend +/// @param runRange Run range +void updateRange(std::pair& selectedTSRange, const std::pair qcTSRange, const std::pair runRange) +{ + if (selectedTSRange.first == qcTSRange.first) { + selectedTSRange.first = runRange.first; + } + if (selectedTSRange.second == qcTSRange.second) { + selectedTSRange.second = runRange.second; + } +} + /// @brief Find bad channels from the occupancy histograms /// @param hits Occupancy histogram /// @param infos Mapping @@ -186,72 +233,61 @@ std::vector getRejectList(std::vector return badChannels; } -/// @brief Gets the run duration with a safety marging -/// @param ccdbApi CCDB api -/// @param marging margin in milliseconds -/// @return Pair with the timestamps of start-margin and end+margin for the run -std::pair getRunDuration(const o2::ccdb::CcdbApi& ccdbApi, int runNumber, int64_t margin = 120000) -{ - auto runRange = o2::ccdb::BasicCCDBManager::getRunDuration(ccdbApi, runNumber); - runRange.first -= margin; // Subtract margin - runRange.second += margin; // Add margin - return runRange; -} - /// @brief Builds the reject list for the selected timestamp -/// @param timestamp Timestamp for query +/// @param md MD structure /// @param qcdbApi QCDB api /// @param ccdbApi CCDB api /// @param outCCDBApi api of the CCDB where the reject list will be uploaded /// @return Reject list -RejectListStruct build_rejectlist(long timestamp, const o2::ccdb::CcdbApi& qcdbApi, const o2::ccdb::CcdbApi& ccdbApi) +RejectListStruct build_rejectlist(const MDStruct& md, const o2::ccdb::CcdbApi& qcdbApi, const o2::ccdb::CcdbApi& ccdbApi) { - std::map metadata; RejectListStruct rl; - auto* qcQuality = qcdbApi.retrieveFromTFileAny(sPathQCQuality, metadata, getTSMS(timestamp)); + if (md.runType != "PHYSICS") { + std::cout << "Run " << md.runNumber << " is of type " << md.runType << ": skip" << std::endl; + return rl; + } + + std::map metadata; + auto* qcQuality = qcdbApi.retrieveFromTFileAny(sPathQCQuality, metadata, getTSMS(md.start)); if (!qcQuality) { - std::cerr << "Cannot find QC quality for " << tsToString(timestamp) << std::endl; + std::cerr << "Cannot find QC quality for " << tsToString(md.start) << std::endl; return rl; } + // Find the first and last timestamp where the quality was bad (if any) auto badTSRange = findTSRange(qcQuality); if (badTSRange.second == 0) { std::cout << "All good" << std::endl; return rl; } + + // Find the first and last timestamp where the quality flag was set + auto qualityTSRange = getFirstLast(qcQuality); // Search for the last timestamp for which the run quality was good auto goodTSRange = findTSRange(qcQuality, false); - // Query the CCDB to see to which run the timestamp corresponds - auto oldestTSInQCQuality = (goodTSRange.first == 0) ? badTSRange.first : goodTSRange.first; - auto grpecs = *ccdbApi.retrieveFromTFileAny("GLO/Config/GRPECS", metadata, getTSMS(oldestTSInQCQuality)); - if (!grpecs.isDetReadOut(o2::detectors::DetID::MID)) { - std::cout << "Error: we are probably reading a parallel run" << std::endl; - grpecs.print(); - return rl; - } - if (grpecs.getRunType() != o2::parameters::GRPECS::PHYSICS) { - std::cout << "This is not a physics run: skip" << std::endl; - grpecs.print(); - return rl; - } - auto runRange = getRunDuration(ccdbApi, grpecs.getRun()); + auto runRange = o2::ccdb::BasicCCDBManager::getRunDuration(ccdbApi, md.runNumber); + updateRange(badTSRange, qualityTSRange, runRange); + updateRange(goodTSRange, qualityTSRange, runRange); // Search for hits histogram in the period where the QC quality was bad - auto tsVector = findObjectsTSInPeriod(badTSRange.first, badTSRange.second, qcdbApi, "qc/MID/MO/QcTaskMIDDigits/Hits"); - if (tsVector.empty()) { + auto mdVector = findObjectsMDInPeriod(badTSRange.first, badTSRange.second, qcdbApi, "qc/MID/MO/QcTaskMIDDigits/Hits"); + if (mdVector.empty()) { std::cerr << "Cannot find hits in period " << tsToString(badTSRange.first) << " - " << tsToString(badTSRange.second) << std::endl; return {}; } - // Focus on the first object found - TH1* occupancy = qcdbApi.retrieveFromTFileAny("qc/MID/MO/QcTaskMIDDigits/Hits", metadata, getTSMS(tsVector.front())); + // Focus on the last object found + // We chose the last instead of the first because it might happen that + // we lose additional boards before the EOR + // If we build the reject list for the first object, we would therefore miss some boards + TH1* occupancy = qcdbApi.retrieveFromTFileAny("qc/MID/MO/QcTaskMIDDigits/Hits", metadata, getTSMS(mdVector.back().start)); o2::mid::GlobalMapper gm; auto infos = gm.buildStripsInfo(); auto badChannels = findBadChannels(occupancy, infos); - auto badChannelsCCDB = *ccdbApi.retrieveFromTFileAny>("MID/Calib/BadChannels", metadata, getTSMS(timestamp)); + auto badChannelsCCDB = *ccdbApi.retrieveFromTFileAny>("MID/Calib/BadChannels", metadata, getTSMS(md.start)); rl.rejectList = getRejectList(badChannels, badChannelsCCDB); if (rl.rejectList.empty()) { - std::cout << "Warning: reject list was empty. It probably means that an entire board is already masked in calibration for run " << grpecs.getRun() << std::endl; + std::cout << "Warning: reject list was empty. It probably means that an entire board is already masked in calibration for run " << md.runNumber << std::endl; return rl; } @@ -260,21 +296,15 @@ RejectListStruct build_rejectlist(long timestamp, const o2::ccdb::CcdbApi& qcdbA for (auto& col : rl.rejectList) { std::cout << col << std::endl; } - std::cout << "Run number: " << grpecs.getRun() << std::endl; - std::cout << "SOR - EOR: " << timeRangeToString(grpecs.getTimeStart(), grpecs.getTimeEnd()) << std::endl; + std::cout << "Run number: " << md.runNumber << std::endl; std::cout << "SOT - EOT: " << timeRangeToString(runRange.first, runRange.second) << std::endl; std::cout << "Good: " << timeRangeToString(goodTSRange.first, goodTSRange.second) << std::endl; std::cout << "Bad: " << timeRangeToString(badTSRange.first, badTSRange.second) << std::endl; + std::cout << "Fraction bad: " << static_cast(badTSRange.second - badTSRange.first) / static_cast(runRange.second - runRange.first) << std::endl; // Set the start of the reject list to the last timestamp in which the occupancy was ok rl.start = goodTSRange.second; - if (goodTSRange.first == 0) { - // If the quality was bad for the full run, set the start of the reject list to the SOR - std::cout << "CAVEAT: no good TS found. Will use SOT instead" << std::endl; - rl.start = runRange.first; - } - // Set the end of the reject list to the end of run - rl.end = runRange.second; + rl.end = badTSRange.second; return rl; } @@ -301,8 +331,8 @@ RejectListStruct load_from_json(const o2::ccdb::CcdbApi& ccdbApi, const char* fi std::cerr << "Problem parsing " << filename << std::endl; return rl; } - auto startRange = getRunDuration(ccdbApi, doc["startRun"].GetInt()); - auto endRange = getRunDuration(ccdbApi, doc["endRun"].GetInt()); + auto startRange = o2::ccdb::BasicCCDBManager::getRunDuration(ccdbApi, doc["startRun"].GetInt()); + auto endRange = o2::ccdb::BasicCCDBManager::getRunDuration(ccdbApi, doc["endRun"].GetInt()); rl.start = startRange.first; rl.end = endRange.second; std::cout << "Manual RL validity: " << timeRangeToString(rl.start, rl.end) << std::endl; @@ -397,9 +427,9 @@ void build_rejectlist(long start, long end, const char* qcdbUrl = "http://ali-qc o2::ccdb::CcdbApi ccdbApi; ccdbApi.init(ccdbUrl); std::vector rls; - auto objectsTS = findObjectsTSInPeriod(start, end, qcdbApi, sPathQCQuality.c_str()); - for (auto ts : objectsTS) { - auto rl = build_rejectlist(ts, qcdbApi, ccdbApi); + auto objectsMD = findObjectsMDInPeriod(start, end, qcdbApi, sPathQCQuality.c_str()); + for (auto md : objectsMD) { + auto rl = build_rejectlist(md, qcdbApi, ccdbApi); if (rl.start != rl.end) { rls.emplace_back(rl); } diff --git a/Detectors/TPC/base/include/TPCBase/CDBTypes.h b/Detectors/TPC/base/include/TPCBase/CDBTypes.h index 27cc2e5a79589..75278f2a76902 100644 --- a/Detectors/TPC/base/include/TPCBase/CDBTypes.h +++ b/Detectors/TPC/base/include/TPCBase/CDBTypes.h @@ -24,68 +24,70 @@ namespace o2::tpc /// Calibration and parameter types for CCDB enum class CDBType { - CalPedestal, ///< Pedestal calibration - CalNoise, ///< Noise calibration - CalPedestalNoise, ///< Pedestal and Noise calibration - CalPulser, ///< Pulser calibration - CalCE, ///< Laser CE calibration - CalPadGainFull, ///< Full pad gain calibration - CalPadGainResidual, ///< ResidualpPad gain calibration (e.g. from tracks) - CalLaserTracks, ///< Laser track calibration data - CalVDriftTgl, ///< ITS-TPC difTgl vdrift calibration - CalTimeGain, ///< Gain variation over time - CalTimeGainMC, ///< Gain variation over time for MC - CalGas, ///< DCS gas measurements - CalTemperature, ///< DCS temperature measurements - CalHV, ///< DCS HV measurements - CalTopologyGain, ///< Q cluster topology correction - /// - ConfigFEEPad, ///< FEE pad-by-pad configuration map - ConfigFEE, ///< FEE configuration map for each tag - ConfigRunInfo, ///< FEE run information (run -> tag) - /// - ParDetector, ///< Parameter for Detector - ParElectronics, ///< Parameter for Electronics - ParGas, ///< Parameter for Gas - ParGEM, ///< Parameter for GEM - /// - CalIDC0A, ///< I_0(r,\phi) = _t - CalIDC0C, ///< I_0(r,\phi) = _t - CalIDC1A, ///< I_1(t) = _{r,\phi} - CalIDC1C, ///< I_1(t) = _{r,\phi} - CalIDCDeltaA, ///< \Delta I(r,\phi,t) = I(r,\phi,t) / ( I_0(r,\phi) * I_1(t) ) - CalIDCDeltaC, ///< \Delta I(r,\phi,t) = I(r,\phi,t) / ( I_0(r,\phi) * I_1(t) ) - CalIDCFourierA, ///< Fourier coefficients of CalIDC1 - CalIDCFourierC, ///< Fourier coefficients of CalIDC1 - CalIDCPadStatusMapA, ///< Status map of the pads (dead etc. obatined from CalIDC0) - CalIDCPadStatusMapC, ///< Status map of the pads (dead etc. obatined from CalIDC0) - CalIDCGroupingParA, ///< Parameters which were used for the averaging of the CalIDCDelta - CalIDCGroupingParC, ///< Parameters which were used for the averaging of the CalIDCDelta - /// - CalSAC0, ///< I_0(r,\phi) = _t - CalSAC1, ///< I_1(t) = _{r,\phi} - CalSACDelta, ///< \Delta I(r,\phi,t) = I(r,\phi,t) / ( I_0(r,\phi) * I_1(t) ) - CalSACFourier, ///< Fourier coefficients of CalSAC1 - /// - CalITPC0, ///< 2D average TPC clusters for longer time interval - CalITPC1, ///< 1D integrated TPC clusters - /// - CalCorrMap, ///< Cluster correction map (high IR rate distortions) - CalCorrMapRef, ///< Cluster correction reference map (static distortions) - CalCorrMapMC, ///< Cluster correction map (high IR rate distortions) for MC - CalCorrDerivMapMC, ///< Cluster correction reference map (static distortions) for MC - /// - CalCorrDerivMap, ///< Cluster correction map (derivative map) - /// - CalTimeSeries, ///< integrated DCAs for longer time interval - CalScaler, ///< Scaler from IDCs or combined estimator - CalScalerWeights, ///< Weights for scalers - CalMShape, ///< calibration object for M-shape distortions - /// - CorrMapParam, ///< parameters for CorrectionMapsLoader configuration - /// - DistortionMapMC, ///< full distortions (static + IR dependant) for MC used in the digitizer - DistortionMapDerivMC ///< derivative distortions for MC used in the digitizer for scaling + CalPedestal, ///< Pedestal calibration + CalNoise, ///< Noise calibration + CalPedestalNoise, ///< Pedestal and Noise calibration + CalPulser, ///< Pulser calibration + CalCE, ///< Laser CE calibration + CalPadGainFull, ///< Full pad gain calibration + CalPadGainResidual, ///< ResidualpPad gain calibration (e.g. from tracks) + CalLaserTracks, ///< Laser track calibration data + CalVDriftTgl, ///< ITS-TPC difTgl vdrift calibration + CalTimeGain, ///< Gain variation over time + CalTimeGainMC, ///< Gain variation over time for MC + CalGas, ///< DCS gas measurements + CalTemperature, ///< DCS temperature measurements + CalHV, ///< DCS HV measurements + CalTopologyGain, ///< Q cluster topology correction + /// + ConfigFEEPad, ///< FEE pad-by-pad configuration map + ConfigFEE, ///< FEE configuration map for each tag + ConfigRunInfo, ///< FEE run information (run -> tag) + /// + ParDetector, ///< Parameter for Detector + ParElectronics, ///< Parameter for Electronics + ParGas, ///< Parameter for Gas + ParGEM, ///< Parameter for GEM + /// + CalIDC0A, ///< I_0(r,\phi) = _t + CalIDC0C, ///< I_0(r,\phi) = _t + CalIDC1A, ///< I_1(t) = _{r,\phi} + CalIDC1C, ///< I_1(t) = _{r,\phi} + CalIDCDeltaA, ///< \Delta I(r,\phi,t) = I(r,\phi,t) / ( I_0(r,\phi) * I_1(t) ) + CalIDCDeltaC, ///< \Delta I(r,\phi,t) = I(r,\phi,t) / ( I_0(r,\phi) * I_1(t) ) + CalIDCFourierA, ///< Fourier coefficients of CalIDC1 + CalIDCFourierC, ///< Fourier coefficients of CalIDC1 + CalIDCPadStatusMapA, ///< Status map of the pads (dead etc. obatined from CalIDC0) + CalIDCPadStatusMapC, ///< Status map of the pads (dead etc. obatined from CalIDC0) + CalIDCGroupingParA, ///< Parameters which were used for the averaging of the CalIDCDelta + CalIDCGroupingParC, ///< Parameters which were used for the averaging of the CalIDCDelta + /// + CalSAC0, ///< I_0(r,\phi) = _t + CalSAC1, ///< I_1(t) = _{r,\phi} + CalSACDelta, ///< \Delta I(r,\phi,t) = I(r,\phi,t) / ( I_0(r,\phi) * I_1(t) ) + CalSACFourier, ///< Fourier coefficients of CalSAC1 + /// + CalITPC0, ///< 2D average TPC clusters for longer time interval + CalITPC1, ///< 1D integrated TPC clusters + /// + CalCorrMap, ///< Cluster correction map (high IR rate distortions) + CalCorrMapRef, ///< Cluster correction reference map (static distortions) + CalCorrMapMC, ///< Cluster correction map (high IR rate distortions) for MC + CalCorrDerivMapMC, ///< Cluster correction reference map (static distortions) for MC + /// + CalCorrDerivMap, ///< Cluster correction map (derivative map) + /// + CalTimeSeries, ///< integrated DCAs for longer time interval + CalScaler, ///< Scaler from IDCs or combined estimator + CalScalerWeights, ///< Weights for scalers + CalMShape, ///< calibration object for M-shape distortions + /// + CorrMapParam, ///< parameters for CorrectionMapsLoader configuration + /// + DistortionMapMC, ///< full distortions (static + IR dependant) for MC used in the digitizer + DistortionMapDerivMC, ///< derivative distortions for MC used in the digitizer for scaling + + AltroSyncSignal ///< timing of Altro chip sync. signal }; /// Storage name in CCDB for each calibration and parameter type @@ -153,6 +155,8 @@ const std::unordered_map CDBTypeMap{ // distortion maps {CDBType::DistortionMapMC, "TPC/Calib/DistortionMapMC"}, {CDBType::DistortionMapDerivMC, "TPC/Calib/DistortionMapDerivativeMC"}, + // AltroSyncSignal + {CDBType::AltroSyncSignal, "TPC/Config/AltroSyncSignal"}, }; } // namespace o2::tpc diff --git a/Detectors/TPC/simulation/src/Detector.cxx b/Detectors/TPC/simulation/src/Detector.cxx index e261424c41332..1a7c0fc25802b 100644 --- a/Detectors/TPC/simulation/src/Detector.cxx +++ b/Detectors/TPC/simulation/src/Detector.cxx @@ -142,6 +142,25 @@ Bool_t Detector::ProcessHits(FairVolume* vol) // TODO: Temporary hack to process only one sector // if (sectorID != 0) return kFALSE; + // ---| momentum and beta gamma |--- + static TLorentzVector momentum; // static to make avoid creation/deletion of this expensive object + fMC->TrackMomentum(momentum); + + const float time = fMC->TrackTime() * 1.0e9; + const int trackID = fMC->GetStack()->GetCurrentTrackNumber(); + const int detID = vol->getMCid(); + o2::data::Stack* stack = (o2::data::Stack*)fMC->GetStack(); + if (fMC->IsTrackEntering() || fMC->IsTrackExiting()) { + stack->addTrackReference(o2::TrackReference(position.X(), position.Y(), position.Z(), momentum.X(), momentum.Y(), + momentum.Z(), fMC->TrackLength(), time, trackID, GetDetId())); + lastReferenceR = fMC->TrackLength(); + } + if (TMath::Abs(lastReferenceR - fMC->TrackLength()) > kMaxDistRef) { /// we can speedup + stack->addTrackReference(o2::TrackReference(position.X(), position.Y(), position.Z(), momentum.X(), momentum.Y(), + momentum.Z(), fMC->TrackLength(), time, trackID, GetDetId())); + lastReferenceR = fMC->TrackLength(); + } + // ---| remove clusters between the IFC and the FC strips |--- // those should not enter the active readout area // do coarse selection before, to limit number of transformations @@ -164,24 +183,6 @@ Bool_t Detector::ProcessHits(FairVolume* vol) } } - // ---| momentum and beta gamma |--- - static TLorentzVector momentum; // static to make avoid creation/deletion of this expensive object - fMC->TrackMomentum(momentum); - - const float time = fMC->TrackTime() * 1.0e9; - const int trackID = fMC->GetStack()->GetCurrentTrackNumber(); - const int detID = vol->getMCid(); - o2::data::Stack* stack = (o2::data::Stack*)fMC->GetStack(); - if (fMC->IsTrackEntering() || fMC->IsTrackExiting()) { - stack->addTrackReference(o2::TrackReference(position.X(), position.Y(), position.Z(), momentum.X(), momentum.Y(), - momentum.Z(), fMC->TrackLength(), time, trackID, GetDetId())); - } - if (TMath::Abs(lastReferenceR - fMC->TrackLength()) > kMaxDistRef) { /// we can speedup - stack->addTrackReference(o2::TrackReference(position.X(), position.Y(), position.Z(), momentum.X(), momentum.Y(), - momentum.Z(), fMC->TrackLength(), time, trackID, GetDetId())); - lastReferenceR = fMC->TrackLength(); - } - // ===| CONVERT THE ENERGY LOSS TO IONIZATION ELECTRONS |===================== // // The energy loss is implemented directly below and taken GEANT3, diff --git a/Framework/Core/include/Framework/AnalysisDataModel.h b/Framework/Core/include/Framework/AnalysisDataModel.h index c90e46bf6da06..e277925ed5603 100644 --- a/Framework/Core/include/Framework/AnalysisDataModel.h +++ b/Framework/Core/include/Framework/AnalysisDataModel.h @@ -15,6 +15,7 @@ #include #include +#include #include #include // std::move @@ -667,28 +668,54 @@ using FullTrack = FullTracks::iterator; namespace trackqa { // TRACKQA TABLE COLUMNS -DECLARE_SOA_INDEX_COLUMN(Track, track); //! track to which this QA information belongs -DECLARE_SOA_COLUMN(TPCTime0, tpcTime0, float); //! tpc only time0 (mTime0 in TPC track) -DECLARE_SOA_COLUMN(TPCDCAR, tpcdcaR, int16_t); //! tpc only DCAr -DECLARE_SOA_COLUMN(TPCDCAZ, tpcdcaZ, int16_t); //! tpc only DCAz -DECLARE_SOA_COLUMN(TPCClusterByteMask, tpcClusterByteMask, uint8_t); //! tracklet bitmask - track defining 8 tracklets (152=8*19 rows) bit set if nCluster>thr (default 5) -DECLARE_SOA_COLUMN(TPCdEdxMax0R, tpcdEdxMax0R, uint8_t); //! TPC dEdxQMax -ROC0/dEdx -DECLARE_SOA_COLUMN(TPCdEdxMax1R, tpcdEdxMax1R, uint8_t); //! TPC dEdxQMax -ROC1/dEdx -DECLARE_SOA_COLUMN(TPCdEdxMax2R, tpcdEdxMax2R, uint8_t); //! TPC dEdxQMax -ROC2/dEdx -DECLARE_SOA_COLUMN(TPCdEdxMax3R, tpcdEdxMax3R, uint8_t); //! TPC dEdxQMax -ROC3/dEdx -DECLARE_SOA_COLUMN(TPCdEdxTot0R, tpcdEdxTot0R, uint8_t); //! TPC dEdxQtot -ROC0/dEdx -DECLARE_SOA_COLUMN(TPCdEdxTot1R, tpcdEdxTot1R, uint8_t); //! TPC dEdxQtot -ROC1/dEdx -DECLARE_SOA_COLUMN(TPCdEdxTot2R, tpcdEdxTot2R, uint8_t); //! TPC dEdxQtot -ROC2/dEdx -DECLARE_SOA_COLUMN(TPCdEdxTot3R, tpcdEdxTot3R, uint8_t); //! TPC dEdxQtot -ROC3/dEdx +DECLARE_SOA_INDEX_COLUMN(Track, track); //! track to which this QA information belongs +DECLARE_SOA_COLUMN(TPCTime0, tpcTime0, float); //! tpc only time0 (mTime0 in TPC track) +DECLARE_SOA_COLUMN(TPCDCAR, tpcdcaR, int16_t); //! tpc only DCAr +DECLARE_SOA_COLUMN(TPCDCAZ, tpcdcaZ, int16_t); //! tpc only DCAz +DECLARE_SOA_COLUMN(TPCClusterByteMask, tpcClusterByteMask, uint8_t); //! tracklet bitmask - track defining 8 tracklets (152=8*19 rows) bit set if nCluster>thr (default 5) +DECLARE_SOA_COLUMN(TPCdEdxMax0R, tpcdEdxMax0R, uint8_t); //! TPC dEdxQMax -ROC0/dEdx +DECLARE_SOA_COLUMN(TPCdEdxMax1R, tpcdEdxMax1R, uint8_t); //! TPC dEdxQMax -ROC1/dEdx +DECLARE_SOA_COLUMN(TPCdEdxMax2R, tpcdEdxMax2R, uint8_t); //! TPC dEdxQMax -ROC2/dEdx +DECLARE_SOA_COLUMN(TPCdEdxMax3R, tpcdEdxMax3R, uint8_t); //! TPC dEdxQMax -ROC3/dEdx +DECLARE_SOA_COLUMN(TPCdEdxTot0R, tpcdEdxTot0R, uint8_t); //! TPC dEdxQtot -ROC0/dEdx +DECLARE_SOA_COLUMN(TPCdEdxTot1R, tpcdEdxTot1R, uint8_t); //! TPC dEdxQtot -ROC1/dEdx +DECLARE_SOA_COLUMN(TPCdEdxTot2R, tpcdEdxTot2R, uint8_t); //! TPC dEdxQtot -ROC2/dEdx +DECLARE_SOA_COLUMN(TPCdEdxTot3R, tpcdEdxTot3R, uint8_t); //! TPC dEdxQtot -ROC3/dEdx +DECLARE_SOA_COLUMN(DeltaRefContParamY, deltaRefContParamY, int8_t); //! Normalized delta of contributor tracks at reference point in the same frame Y +DECLARE_SOA_COLUMN(DeltaRefContParamZ, deltaRefITSParamZ, int8_t); //! Normalized delta of contributor tracks at reference point in the same frame Z +DECLARE_SOA_COLUMN(DeltaRefContParamSnp, deltaRefContParamSnp, int8_t); //! Normalized delta of contributor tracks at reference point in the same frame Snp +DECLARE_SOA_COLUMN(DeltaRefContParamTgl, deltaRefContParamTgl, int8_t); //! Normalized delta of contributor tracks at reference point in the same frame Tgl +DECLARE_SOA_COLUMN(DeltaRefContParamQ2Pt, deltaRefContParamQ2Pt, int8_t); //! Normalized delta of contributor tracks at reference point in the same frame Q2Pt +DECLARE_SOA_COLUMN(DeltaRefGloParamY, deltaRefGloParamY, int8_t); //! Normalized delta of global track to average contributors matched tracks at reference point in the same frame Y +DECLARE_SOA_COLUMN(DeltaRefGloParamZ, deltaRefGloParamZ, int8_t); //! Normalized delta of global track to average contributors matched tracks at reference point in the same frame Z +DECLARE_SOA_COLUMN(DeltaRefGloParamSnp, deltaRefGloParamSnp, int8_t); //! Normalized delta of global track to average contributors matched tracks at reference point in the same frame Snp +DECLARE_SOA_COLUMN(DeltaRefGloParamTgl, deltaRefGloParamTgl, int8_t); //! Normalized delta of global track to average contributors matched tracks at reference point in the same frame Tgl +DECLARE_SOA_COLUMN(DeltaRefGloParamQ2Pt, deltaRefGloParamQ2Pt, int8_t); //! Normalized delta of global track to average contributors matched tracks at reference point in the same frame Q2Pt + +DECLARE_SOA_DYNAMIC_COLUMN(IsDummy, isDummy, //! indicates if the propagation of the contrib. tracks was successful and residuals are available + [](int8_t cY, int8_t cZ, int8_t cSnp, int8_t cTgl, int8_t cQ2Pt, int8_t gY, int8_t gZ, int8_t gSnp, int8_t gTgl, int8_t gQ2Pt) -> bool { + constexpr int8_t m = std::numeric_limits::min(); + return (cY == m && cZ == m && cSnp == m && cTgl == m && cQ2Pt == m && gY == m && gZ == m && gSnp == m && gTgl == m && gQ2Pt == m); + }); } // namespace trackqa -DECLARE_SOA_TABLE(TracksQA, "AOD", "TRACKQA", //! trackQA information - sampled QA information currently for the TPC +DECLARE_SOA_TABLE(TracksQA_000, "AOD", "TRACKQA", //! trackQA information - sampled QA information currently for the TPC - version 0 o2::soa::Index<>, trackqa::TrackId, trackqa::TPCTime0, trackqa::TPCDCAR, trackqa::TPCDCAZ, trackqa::TPCClusterByteMask, trackqa::TPCdEdxMax0R, trackqa::TPCdEdxMax1R, trackqa::TPCdEdxMax2R, trackqa::TPCdEdxMax3R, trackqa::TPCdEdxTot0R, trackqa::TPCdEdxTot1R, trackqa::TPCdEdxTot2R, trackqa::TPCdEdxTot3R); // o2::soa::Index<>, trackqa::TrackId, trackqa::TPCDCAR, trackqa::TPCDCAZ, trackqa::TPCClusterByteMask, -using TrackQA = TracksQA::iterator; +DECLARE_SOA_TABLE_VERSIONED(TracksQA_001, "AOD", "TRACKQA", 1, //! trackQA information - version 1 - including contributor residuals of matched tracks at reference radius + o2::soa::Index<>, trackqa::TrackId, trackqa::TPCTime0, trackqa::TPCDCAR, trackqa::TPCDCAZ, trackqa::TPCClusterByteMask, + trackqa::TPCdEdxMax0R, trackqa::TPCdEdxMax1R, trackqa::TPCdEdxMax2R, trackqa::TPCdEdxMax3R, + trackqa::TPCdEdxTot0R, trackqa::TPCdEdxTot1R, trackqa::TPCdEdxTot2R, trackqa::TPCdEdxTot3R, + trackqa::DeltaRefContParamY, trackqa::DeltaRefContParamZ, trackqa::DeltaRefContParamSnp, trackqa::DeltaRefContParamTgl, trackqa::DeltaRefContParamQ2Pt, + trackqa::DeltaRefGloParamY, trackqa::DeltaRefGloParamZ, trackqa::DeltaRefGloParamSnp, trackqa::DeltaRefGloParamTgl, trackqa::DeltaRefGloParamQ2Pt, + trackqa::IsDummy); + +using TracksQAVersion = TracksQA_000; +using TracksQA = TracksQAVersion::iterator; namespace fwdtrack { diff --git a/Framework/Core/include/Framework/DataTypes.h b/Framework/Core/include/Framework/DataTypes.h index 92af1f79e2314..9d829159718d8 100644 --- a/Framework/Core/include/Framework/DataTypes.h +++ b/Framework/Core/include/Framework/DataTypes.h @@ -15,6 +15,7 @@ #include #include +#include namespace o2::aod::bc { @@ -120,6 +121,15 @@ struct TPCTimeErrEncoding { } }; } // namespace extensions + +// Reference radius for extrapolated tracks +constexpr float trackQARefRadius{50.f}; +constexpr float trackQAScaleBins{5.f}; +// Fit parameters for scale dY, dZ, dSnp, dTgl, dQ2Pt +constexpr std::array trackQAScaleContP0{0.257192, 0.0775375, 0.00424283, 0.00107201, 0.0335447}; +constexpr std::array trackQAScaleContP1{0.189371, 0.409071, 0.00694444, 0.00720038, 0.0806902}; +constexpr std::array trackQAScaleGloP0{0.130985, 0.0775375, 0.00194703, 0.000405458, 0.0160007}; +constexpr std::array trackQAScaleGloP1{0.183731, 0.409071, 0.00621802, 0.00624881, 0.0418957}; } // namespace o2::aod::track namespace o2::aod::fwdtrack diff --git a/Framework/Core/src/DataDescriptorQueryBuilder.cxx b/Framework/Core/src/DataDescriptorQueryBuilder.cxx index 41a14d06f3acc..8b0c239699cc9 100644 --- a/Framework/Core/src/DataDescriptorQueryBuilder.cxx +++ b/Framework/Core/src/DataDescriptorQueryBuilder.cxx @@ -319,12 +319,9 @@ std::vector DataDescriptorQueryBuilder::parse(char const* config) if (*currentKey == "lifetime" && currentValue == "condition") { currentLifetime = Lifetime::Condition; } - if (*currentKey == "ccdb-run-dependent" && (currentValue != "false" && currentValue != "0")) { - attributes.push_back(ConfigParamSpec{*currentKey, VariantType::Bool, true, {}}); - } else if (*currentKey == "ccdb-run-dependent" && (currentValue == "false" || currentValue == "0")) { - attributes.push_back(ConfigParamSpec{*currentKey, VariantType::Bool, false, {}}); - } else if (*currentKey == "ccdb-run-dependent") { - error("ccdb-run-dependent can only be true or false"); + if (*currentKey == "ccdb-run-dependent") { + int val = currentValue == "false" ? 0 : (currentValue == "true" ? 1 : std::stoi(*currentValue)); + attributes.push_back(ConfigParamSpec{*currentKey, VariantType::Int, val, {}}); } else { attributes.push_back(ConfigParamSpec{*currentKey, VariantType::String, *currentValue, {}}); } @@ -333,12 +330,9 @@ std::vector DataDescriptorQueryBuilder::parse(char const* config) if (*currentKey == "lifetime" && currentValue == "condition") { currentLifetime = Lifetime::Condition; } - if (*currentKey == "ccdb-run-dependent" && (currentValue != "false" && currentValue != "0")) { - attributes.push_back(ConfigParamSpec{*currentKey, VariantType::Bool, true, {}}); - } else if (*currentKey == "ccdb-run-dependent" && (currentValue == "false" || currentValue == "0")) { - attributes.push_back(ConfigParamSpec{*currentKey, VariantType::Bool, false, {}}); - } else if (*currentKey == "ccdb-run-dependent") { - error("ccdb-run-dependent can only be true or false"); + if (*currentKey == "ccdb-run-dependent") { + int val = currentValue == "false" ? 0 : (currentValue == "true" ? 1 : std::stoi(*currentValue)); + attributes.push_back(ConfigParamSpec{*currentKey, VariantType::Int, val, {}}); } else { attributes.push_back(ConfigParamSpec{*currentKey, VariantType::String, *currentValue, {}}); } @@ -347,12 +341,9 @@ std::vector DataDescriptorQueryBuilder::parse(char const* config) if (*currentKey == "lifetime" && currentValue == "condition") { currentLifetime = Lifetime::Condition; } - if (*currentKey == "ccdb-run-dependent" && (currentValue != "false" && currentValue != "0")) { - attributes.push_back(ConfigParamSpec{*currentKey, VariantType::Bool, true, {}}); - } else if (*currentKey == "ccdb-run-dependent" && (currentValue == "false" || currentValue == "0")) { - attributes.push_back(ConfigParamSpec{*currentKey, VariantType::Bool, false, {}}); - } else if (*currentKey == "ccdb-run-dependent") { - error("ccdb-run-dependent can only be true or false"); + if (*currentKey == "ccdb-run-dependent") { + int val = currentValue == "false" ? 0 : (currentValue == "true" ? 1 : std::stoi(*currentValue)); + attributes.push_back(ConfigParamSpec{*currentKey, VariantType::Int, val, {}}); } else { attributes.push_back(ConfigParamSpec{*currentKey, VariantType::String, *currentValue, {}}); } diff --git a/Framework/Core/src/RootArrowFilesystem.cxx b/Framework/Core/src/RootArrowFilesystem.cxx index 7e331814272a6..5f2d21d942d37 100644 --- a/Framework/Core/src/RootArrowFilesystem.cxx +++ b/Framework/Core/src/RootArrowFilesystem.cxx @@ -13,6 +13,7 @@ #include "Framework/RuntimeError.h" #include "Framework/Signpost.h" #include +#include #include #include #include @@ -427,7 +428,7 @@ class TTreeFileWriter : public arrow::dataset::FileWriter O2_SIGNPOST_EVENT_EMIT(root_arrow_fs, sid, "finaliseBasketSize", "Branch %s exists and uses %d bytes per entry.", branch->GetName(), valueSize); // This should probably lookup the - auto column = firstBatch->GetColumnByName(branch->GetName()); + auto column = firstBatch->GetColumnByName(schema_->field(i)->name()); auto list = std::static_pointer_cast(column); O2_SIGNPOST_EVENT_EMIT(root_arrow_fs, sid, "finaliseBasketSize", "Branch %s needed. Associated size branch %s and there are %lli entries of size %d in that list.", branch->GetName(), sizeBranch->GetName(), list->length(), valueSize); @@ -497,8 +498,8 @@ class TTreeFileWriter : public arrow::dataset::FileWriter } break; case arrow::Type::LIST: { valueTypes.push_back(field->type()->field(0)->type()); - listSizes.back() = 0; // VLA, we need to calculate it on the fly; std::string leafList = fmt::format("{}[{}_size]{}", field->name(), field->name(), rootSuffixFromArrow(valueTypes.back()->id())); + listSizes.back() = -1; // VLA, we need to calculate it on the fly; std::string sizeLeafList = field->name() + "_size/I"; sizesBranches.push_back(treeStream->CreateBranch((field->name() + "_size").c_str(), sizeLeafList.c_str())); branches.push_back(treeStream->CreateBranch(field->name().c_str(), leafList.c_str())); @@ -765,7 +766,7 @@ arrow::Result TTreeFileFormat::ScanBatchesAsync( typeSize = fixedSizeList->field(0)->type()->byte_width(); } else if (auto vlaListType = std::dynamic_pointer_cast(physicalField->type())) { listSize = -1; - typeSize = fixedSizeList->field(0)->type()->byte_width(); + typeSize = vlaListType->field(0)->type()->byte_width(); } if (listSize == -1) { mSizeBranch = branch->GetTree()->GetBranch((std::string{branch->GetName()} + "_size").c_str()); diff --git a/Framework/Core/test/test_Root2ArrowTable.cxx b/Framework/Core/test/test_Root2ArrowTable.cxx index a659d488ae24a..2b0ab9154250c 100644 --- a/Framework/Core/test/test_Root2ArrowTable.cxx +++ b/Framework/Core/test/test_Root2ArrowTable.cxx @@ -322,12 +322,26 @@ bool validateContents(std::shared_ptr batch) REQUIRE(bool_array->Value(1) == (i % 5 == 0)); } } + + { + auto list_array = std::static_pointer_cast(batch->GetColumnByName("vla")); + + REQUIRE(list_array->length() == 100); + for (int64_t i = 0; i < list_array->length(); i++) { + auto value_slice = list_array->value_slice(i); + REQUIRE(value_slice->length() == (i % 10)); + auto int_array = std::static_pointer_cast(value_slice); + for (size_t j = 0; j < value_slice->length(); j++) { + REQUIRE(int_array->Value(j) == j); + } + } + } return true; } bool validateSchema(std::shared_ptr schema) { - REQUIRE(schema->num_fields() == 9); + REQUIRE(schema->num_fields() == 10); REQUIRE(schema->field(0)->type()->id() == arrow::float32()->id()); REQUIRE(schema->field(1)->type()->id() == arrow::float32()->id()); REQUIRE(schema->field(2)->type()->id() == arrow::float32()->id()); @@ -337,6 +351,7 @@ bool validateSchema(std::shared_ptr schema) REQUIRE(schema->field(6)->type()->id() == arrow::fixed_size_list(arrow::int32(), 2)->id()); REQUIRE(schema->field(7)->type()->id() == arrow::boolean()->id()); REQUIRE(schema->field(8)->type()->id() == arrow::fixed_size_list(arrow::boolean(), 2)->id()); + REQUIRE(schema->field(9)->type()->id() == arrow::list(arrow::int32())->id()); return true; } @@ -390,6 +405,8 @@ TEST_CASE("RootTree2Dataset") Int_t ev; bool oneBool; bool manyBool[2]; + int vla[10] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9}; + int vlaSize = 0; t->Branch("px", &px, "px/F"); t->Branch("py", &py, "py/F"); @@ -400,6 +417,8 @@ TEST_CASE("RootTree2Dataset") t->Branch("ij", ij, "ij[2]/I"); t->Branch("bools", &oneBool, "bools/O"); t->Branch("manyBools", &manyBool, "manyBools[2]/O"); + t->Branch("vla_size", &vlaSize, "vla_size/I"); + t->Branch("vla", vla, "vla[vla_size]/I"); // fill the tree for (Int_t i = 0; i < 100; i++) { xyz[0] = 1; @@ -415,9 +434,11 @@ TEST_CASE("RootTree2Dataset") oneBool = (i % 3 == 0); manyBool[0] = (i % 4 == 0); manyBool[1] = (i % 5 == 0); + vlaSize = i % 10; t->Fill(); } } + f->Write(); size_t totalSizeCompressed = 0; size_t totalSizeUncompressed = 0; @@ -428,16 +449,7 @@ TEST_CASE("RootTree2Dataset") auto schemaOpt = format->Inspect(source); REQUIRE(schemaOpt.ok()); auto schema = *schemaOpt; - REQUIRE(schema->num_fields() == 9); - REQUIRE(schema->field(0)->type()->id() == arrow::float32()->id()); - REQUIRE(schema->field(1)->type()->id() == arrow::float32()->id()); - REQUIRE(schema->field(2)->type()->id() == arrow::float32()->id()); - REQUIRE(schema->field(3)->type()->id() == arrow::float64()->id()); - REQUIRE(schema->field(4)->type()->id() == arrow::int32()->id()); - REQUIRE(schema->field(5)->type()->id() == arrow::fixed_size_list(arrow::float32(), 3)->id()); - REQUIRE(schema->field(6)->type()->id() == arrow::fixed_size_list(arrow::int32(), 2)->id()); - REQUIRE(schema->field(7)->type()->id() == arrow::boolean()->id()); - REQUIRE(schema->field(8)->type()->id() == arrow::fixed_size_list(arrow::boolean(), 2)->id()); + validateSchema(schema); auto fragment = format->MakeFragment(source, {}, schema); REQUIRE(fragment.ok()); @@ -448,41 +460,9 @@ TEST_CASE("RootTree2Dataset") auto batches = (*scanner)(); auto result = batches.result(); REQUIRE(result.ok()); - REQUIRE((*result)->columns().size() == 9); + REQUIRE((*result)->columns().size() == 10); REQUIRE((*result)->num_rows() == 100); - - { - auto int_array = std::static_pointer_cast((*result)->GetColumnByName("ev")); - for (int64_t j = 0; j < int_array->length(); j++) { - REQUIRE(int_array->Value(j) == j + 1); - } - } - - { - auto list_array = std::static_pointer_cast((*result)->GetColumnByName("xyz")); - - // Iterate over the FixedSizeListArray - for (int64_t i = 0; i < list_array->length(); i++) { - auto value_slice = list_array->value_slice(i); - auto float_array = std::static_pointer_cast(value_slice); - - REQUIRE(float_array->Value(0) == 1); - REQUIRE(float_array->Value(1) == 2); - REQUIRE(float_array->Value(2) == i + 1); - } - } - - { - auto list_array = std::static_pointer_cast((*result)->GetColumnByName("ij")); - - // Iterate over the FixedSizeListArray - for (int64_t i = 0; i < list_array->length(); i++) { - auto value_slice = list_array->value_slice(i); - auto int_array = std::static_pointer_cast(value_slice); - REQUIRE(int_array->Value(0) == i); - REQUIRE(int_array->Value(1) == i + 1); - } - } + validateContents(*result); auto* output = new TMemFile("foo", "RECREATE"); auto outFs = std::make_shared(output, 0); @@ -497,31 +477,31 @@ TEST_CASE("RootTree2Dataset") auto success = writer->get()->Write(*result); auto rootDestination = std::dynamic_pointer_cast(*destination); - REQUIRE(success.ok()); - // Let's read it back... - arrow::dataset::FileSource source2("/DF_3", outFs); - auto newTreeFS = outFs->GetSubFilesystem(source2); - - REQUIRE(format->IsSupported(source) == true); - - auto schemaOptWritten = format->Inspect(source); - REQUIRE(schemaOptWritten.ok()); - auto schemaWritten = *schemaOptWritten; - REQUIRE(validateSchema(schemaWritten)); - - auto fragmentWritten = format->MakeFragment(source, {}, schema); - REQUIRE(fragmentWritten.ok()); - auto optionsWritten = std::make_shared(); - options->dataset_schema = schemaWritten; - auto scannerWritten = format->ScanBatchesAsync(optionsWritten, *fragment); - REQUIRE(scannerWritten.ok()); - auto batchesWritten = (*scanner)(); - auto resultWritten = batches.result(); - REQUIRE(resultWritten.ok()); - REQUIRE((*resultWritten)->columns().size() == 9); - REQUIRE((*resultWritten)->num_rows() == 100); - validateContents(*resultWritten); - + SECTION("Read tree") { + REQUIRE(success.ok()); + // Let's read it back... + arrow::dataset::FileSource source2("/DF_3", outFs); + auto newTreeFS = outFs->GetSubFilesystem(source2); + + REQUIRE(format->IsSupported(source) == true); + + auto schemaOptWritten = format->Inspect(source); + REQUIRE(schemaOptWritten.ok()); + auto schemaWritten = *schemaOptWritten; + REQUIRE(validateSchema(schemaWritten)); + + auto fragmentWritten = format->MakeFragment(source, {}, schema); + REQUIRE(fragmentWritten.ok()); + auto optionsWritten = std::make_shared(); + options->dataset_schema = schemaWritten; + auto scannerWritten = format->ScanBatchesAsync(optionsWritten, *fragment); + REQUIRE(scannerWritten.ok()); + auto batchesWritten = (*scanner)(); + auto resultWritten = batches.result(); + REQUIRE(resultWritten.ok()); + REQUIRE((*resultWritten)->columns().size() == 10); + REQUIRE((*resultWritten)->num_rows() == 100); + validateContents(*resultWritten); } } diff --git a/GPU/Common/GPUCommonConstants.h b/GPU/Common/GPUCommonConstants.h index 883f64b7bdd12..f45aa05ed00ca 100644 --- a/GPU/Common/GPUCommonConstants.h +++ b/GPU/Common/GPUCommonConstants.h @@ -20,7 +20,7 @@ #if !defined(__OPENCL1__) namespace GPUCA_NAMESPACE::gpu::gpu_common_constants { -static CONSTEXPR const float kCLight = 0.000299792458f; +static CONSTEXPR const float kCLight = 0.000299792458f; // TODO: Duplicate of MathConstants, fix this when OpenCL1 is removed } #endif diff --git a/GPU/GPUTracking/Base/GPUParam.cxx b/GPU/GPUTracking/Base/GPUParam.cxx index 42d4f61f77116..661ae830ca6f3 100644 --- a/GPU/GPUTracking/Base/GPUParam.cxx +++ b/GPU/GPUTracking/Base/GPUParam.cxx @@ -120,11 +120,12 @@ void GPUParam::SetDefaults(float solenoidBz) par.toyMCEventsFlag = false; par.continuousTracking = false; continuousMaxTimeBin = 0; + tpcCutTimeBin = 0; par.debugLevel = 0; par.earlyTpcTransform = false; } -void GPUParam::UpdateSettings(const GPUSettingsGRP* g, const GPUSettingsProcessing* p, const GPURecoStepConfiguration* w) +void GPUParam::UpdateSettings(const GPUSettingsGRP* g, const GPUSettingsProcessing* p, const GPURecoStepConfiguration* w, const GPUSettingsRecDynamic* d) { if (g) { UpdateBzOnly(g->solenoidBzNominalGPU); @@ -132,6 +133,7 @@ void GPUParam::UpdateSettings(const GPUSettingsGRP* g, const GPUSettingsProcessi par.toyMCEventsFlag = g->homemadeEvents; par.continuousTracking = g->grpContinuousMaxTimeBin != 0; continuousMaxTimeBin = g->grpContinuousMaxTimeBin == -1 ? GPUSettings::TPC_MAX_TF_TIME_BIN : g->grpContinuousMaxTimeBin; + tpcCutTimeBin = g->tpcCutTimeBin; } par.earlyTpcTransform = rec.tpc.forceEarlyTransform == -1 ? (!par.continuousTracking) : rec.tpc.forceEarlyTransform; qptB5Scaler = CAMath::Abs(bzkG) > 0.1f ? CAMath::Abs(bzkG) / 5.006680f : 1.f; // Repeat here, since passing in g is optional @@ -145,6 +147,9 @@ void GPUParam::UpdateSettings(const GPUSettingsGRP* g, const GPUSettingsProcessi dodEdxDownscaled = (rand() % 100) < p->tpcDownscaledEdx; } } + if (d) { + rec.dyn = *d; + } } void GPUParam::UpdateBzOnly(float newSolenoidBz) diff --git a/GPU/GPUTracking/Base/GPUParam.h b/GPU/GPUTracking/Base/GPUParam.h index 070ac76f58ffb..ce9ac30b7c35b 100644 --- a/GPU/GPUTracking/Base/GPUParam.h +++ b/GPU/GPUTracking/Base/GPUParam.h @@ -59,6 +59,7 @@ struct GPUParam_t { int8_t dodEdxDownscaled; int32_t continuousMaxTimeBin; + int32_t tpcCutTimeBin; GPUTPCGeometry tpcGeometry; // TPC Geometry GPUTPCGMPolynomialField polynomialField; // Polynomial approx. of magnetic field for TPC GM @@ -84,7 +85,7 @@ struct GPUParam : public internal::GPUParam_t #ifndef GPUCA_GPUCODE void SetDefaults(float solenoidBz); void SetDefaults(const GPUSettingsGRP* g, const GPUSettingsRec* r = nullptr, const GPUSettingsProcessing* p = nullptr, const GPURecoStepConfiguration* w = nullptr); - void UpdateSettings(const GPUSettingsGRP* g, const GPUSettingsProcessing* p = nullptr, const GPURecoStepConfiguration* w = nullptr); + void UpdateSettings(const GPUSettingsGRP* g, const GPUSettingsProcessing* p = nullptr, const GPURecoStepConfiguration* w = nullptr, const GPUSettingsRecDynamic* d = nullptr); void UpdateBzOnly(float newSolenoidBz); void LoadClusterErrors(bool Print = 0); void UpdateRun3ClusterErrors(const float* yErrorParam, const float* zErrorParam); diff --git a/GPU/GPUTracking/Base/GPUReconstruction.cxx b/GPU/GPUTracking/Base/GPUReconstruction.cxx index 632bf0f331f31..9abe225c7848e 100644 --- a/GPU/GPUTracking/Base/GPUReconstruction.cxx +++ b/GPU/GPUTracking/Base/GPUReconstruction.cxx @@ -427,7 +427,7 @@ int32_t GPUReconstruction::InitPhaseAfterDevice() (mProcessors[i].proc->*(mProcessors[i].InitializeProcessor))(); } - WriteConstantParams(); // First initialization, if the user doesn't use RunChains + WriteConstantParams(); // Initialize with initial values, can optionally be updated later mInitialized = true; return 0; @@ -1105,7 +1105,12 @@ void GPUReconstruction::DumpSettings(const char* dir) } } -void GPUReconstruction::UpdateSettings(const GPUSettingsGRP* g, const GPUSettingsProcessing* p) +void GPUReconstruction::UpdateDynamicSettings(const GPUSettingsRecDynamic* d) +{ + UpdateSettings(nullptr, nullptr, d); +} + +void GPUReconstruction::UpdateSettings(const GPUSettingsGRP* g, const GPUSettingsProcessing* p, const GPUSettingsRecDynamic* d) { if (g) { mGRPSettings = *g; @@ -1114,8 +1119,11 @@ void GPUReconstruction::UpdateSettings(const GPUSettingsGRP* g, const GPUSetting mProcessingSettings.debugLevel = p->debugLevel; mProcessingSettings.resetTimers = p->resetTimers; } - GPURecoStepConfiguration w = mRecoSteps; - param().UpdateSettings(g, p, &w); + GPURecoStepConfiguration* w = nullptr; + if (mRecoSteps.steps.isSet(GPUDataTypes::RecoStep::TPCdEdx)) { + w = &mRecoSteps; + } + param().UpdateSettings(g, p, w, d); if (mInitialized) { WriteConstantParams(); } diff --git a/GPU/GPUTracking/Base/GPUReconstruction.h b/GPU/GPUTracking/Base/GPUReconstruction.h index 70a066532d938..efad0b41fd571 100644 --- a/GPU/GPUTracking/Base/GPUReconstruction.h +++ b/GPU/GPUTracking/Base/GPUReconstruction.h @@ -200,7 +200,8 @@ class GPUReconstruction void SetSettings(const GPUSettingsGRP* grp, const GPUSettingsRec* rec = nullptr, const GPUSettingsProcessing* proc = nullptr, const GPURecoStepConfiguration* workflow = nullptr); void SetResetTimers(bool reset) { mProcessingSettings.resetTimers = reset; } // May update also after Init() void SetDebugLevelTmp(int32_t level) { mProcessingSettings.debugLevel = level; } // Temporarily, before calling SetSettings() - void UpdateSettings(const GPUSettingsGRP* g, const GPUSettingsProcessing* p = nullptr); + void UpdateSettings(const GPUSettingsGRP* g, const GPUSettingsProcessing* p = nullptr, const GPUSettingsRecDynamic* d = nullptr); + void UpdateDynamicSettings(const GPUSettingsRecDynamic* d); void SetOutputControl(const GPUOutputControl& v) { mOutputControl = v; } void SetOutputControl(void* ptr, size_t size); void SetInputControl(void* ptr, size_t size); diff --git a/GPU/GPUTracking/Base/GPUReconstructionCPU.cxx b/GPU/GPUTracking/Base/GPUReconstructionCPU.cxx index 537c3cf63a628..271bee59db31b 100644 --- a/GPU/GPUTracking/Base/GPUReconstructionCPU.cxx +++ b/GPU/GPUTracking/Base/GPUReconstructionCPU.cxx @@ -228,7 +228,7 @@ int32_t GPUReconstructionCPU::RunChains() mThreadId = GetThread(); } if (mSlaves.size() || mMaster) { - WriteConstantParams(); // Reinitialize + WriteConstantParams(); // Reinitialize // TODO: Get this in sync with GPUChainTracking::DoQueuedUpdates, and consider the doublePipeline } for (uint32_t i = 0; i < mChains.size(); i++) { int32_t retVal = mChains[i]->RunChain(); diff --git a/GPU/GPUTracking/DataCompression/GPUTPCDecompression.cxx b/GPU/GPUTracking/DataCompression/GPUTPCDecompression.cxx index 4039ebb0c100d..7c10f0eeef74f 100644 --- a/GPU/GPUTracking/DataCompression/GPUTPCDecompression.cxx +++ b/GPU/GPUTracking/DataCompression/GPUTPCDecompression.cxx @@ -43,7 +43,7 @@ void GPUTPCDecompression::SetPointersCompressedClusters(void*& mem, T& c, uint32 uint32_t nClAreduced = reducedClA ? nClA - nTr : nClA; - if (!(mRec->GetParam().rec.tpc.compressionTypeMask & GPUSettings::CompressionTrackModel)) { + if (!(c.nComppressionModes & GPUSettings::CompressionTrackModel)) { return; // Track model disabled, do not allocate memory } computePointerWithAlignment(mem, c.qTotA, nClA); @@ -84,6 +84,24 @@ void* GPUTPCDecompression::SetPointersTmpNativeBuffersInput(void* mem) return mem; } +void* GPUTPCDecompression::SetPointersTmpClusterNativeAccessForFiltering(void* mem) +{ + computePointerWithAlignment(mem, mNativeClustersBuffer, mNClusterNativeBeforeFiltering); + return mem; +} + +void* GPUTPCDecompression::SetPointersInputClusterNativeAccess(void* mem) +{ + computePointerWithAlignment(mem, mClusterNativeAccess); + return mem; +} + +void* GPUTPCDecompression::SetPointersNClusterPerSectorRow(void* mem) +{ + computePointerWithAlignment(mem, mNClusterPerSectorRow, NSLICES * GPUCA_ROW_COUNT); + return mem; +} + void GPUTPCDecompression::RegisterMemoryAllocation() { AllocateAndInitializeLate(); @@ -91,6 +109,9 @@ void GPUTPCDecompression::RegisterMemoryAllocation() mRec->RegisterMemoryAllocation(this, &GPUTPCDecompression::SetPointersTmpNativeBuffersGPU, GPUMemoryResource::MEMORY_SCRATCH, "TPCDecompressionTmpBuffersGPU"); mResourceTmpIndexes = mRec->RegisterMemoryAllocation(this, &GPUTPCDecompression::SetPointersTmpNativeBuffersOutput, GPUMemoryResource::MEMORY_OUTPUT | GPUMemoryResource::MEMORY_SCRATCH, "TPCDecompressionTmpBuffersOutput"); mResourceTmpClustersOffsets = mRec->RegisterMemoryAllocation(this, &GPUTPCDecompression::SetPointersTmpNativeBuffersInput, GPUMemoryResource::MEMORY_INPUT | GPUMemoryResource::MEMORY_SCRATCH, "TPCDecompressionTmpBuffersInput"); + mResourceTmpBufferBeforeFiltering = mRec->RegisterMemoryAllocation(this, &GPUTPCDecompression::SetPointersTmpClusterNativeAccessForFiltering, GPUMemoryResource::MEMORY_CUSTOM | GPUMemoryResource::MEMORY_SCRATCH, "TPCDecompressionTmpBufferForFiltering"); + mResourceClusterNativeAccess = mRec->RegisterMemoryAllocation(this, &GPUTPCDecompression::SetPointersInputClusterNativeAccess, GPUMemoryResource::MEMORY_INPUT | GPUMemoryResource::MEMORY_CUSTOM | GPUMemoryResource::MEMORY_SCRATCH, "TPCDecompressionTmpClusterAccessForFiltering"); + mResourceNClusterPerSectorRow = mRec->RegisterMemoryAllocation(this, &GPUTPCDecompression::SetPointersNClusterPerSectorRow, GPUMemoryResource::MEMORY_OUTPUT | GPUMemoryResource::MEMORY_CUSTOM | GPUMemoryResource::MEMORY_SCRATCH, "TPCDecompressionTmpClusterCountForFiltering"); } void GPUTPCDecompression::SetMaxData(const GPUTrackingInOutPointers& io) diff --git a/GPU/GPUTracking/DataCompression/GPUTPCDecompression.h b/GPU/GPUTracking/DataCompression/GPUTPCDecompression.h index d9871613d8401..47c64008b176e 100644 --- a/GPU/GPUTracking/DataCompression/GPUTPCDecompression.h +++ b/GPU/GPUTracking/DataCompression/GPUTPCDecompression.h @@ -55,6 +55,9 @@ class GPUTPCDecompression : public GPUProcessor void* SetPointersTmpNativeBuffersGPU(void* mem); void* SetPointersTmpNativeBuffersOutput(void* mem); void* SetPointersTmpNativeBuffersInput(void* mem); + void* SetPointersTmpClusterNativeAccessForFiltering(void* mem); + void* SetPointersInputClusterNativeAccess(void* mem); + void* SetPointersNClusterPerSectorRow(void* mem); #endif @@ -63,11 +66,14 @@ class GPUTPCDecompression : public GPUProcessor o2::tpc::CompressedClusters mInputGPU; uint32_t mMaxNativeClustersPerBuffer; + uint32_t mNClusterNativeBeforeFiltering; uint32_t* mNativeClustersIndex; uint32_t* mUnattachedClustersOffsets; uint32_t* mAttachedClustersOffsets; + uint32_t* mNClusterPerSectorRow; o2::tpc::ClusterNative* mTmpNativeClusters; o2::tpc::ClusterNative* mNativeClustersBuffer; + o2::tpc::ClusterNativeAccess* mClusterNativeAccess; template void SetPointersCompressedClusters(void*& mem, T& c, uint32_t nClA, uint32_t nTr, uint32_t nClU, bool reducedClA); @@ -75,6 +81,9 @@ class GPUTPCDecompression : public GPUProcessor int16_t mMemoryResInputGPU = -1; int16_t mResourceTmpIndexes = -1; int16_t mResourceTmpClustersOffsets = -1; + int16_t mResourceTmpBufferBeforeFiltering = -1; + int16_t mResourceClusterNativeAccess = -1; + int16_t mResourceNClusterPerSectorRow = -1; }; } // namespace GPUCA_NAMESPACE::gpu #endif // GPUTPCDECOMPRESSION_H diff --git a/GPU/GPUTracking/DataCompression/GPUTPCDecompressionKernels.cxx b/GPU/GPUTracking/DataCompression/GPUTPCDecompressionKernels.cxx index 2c88ea0079a26..d7f1e2ac88368 100644 --- a/GPU/GPUTracking/DataCompression/GPUTPCDecompressionKernels.cxx +++ b/GPU/GPUTracking/DataCompression/GPUTPCDecompressionKernels.cxx @@ -43,7 +43,7 @@ GPUdii() void GPUTPCDecompressionKernels::Thread 0 ? cl.getTime() < param.tpcCutTimeBin : true; +} + +template <> +GPUdii() void GPUTPCDecompressionUtilKernels::Thread(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread, GPUsharedref() GPUSharedMemory& smem, processorType& processors) +{ + const GPUParam& GPUrestrict() param = processors.param; + GPUTPCDecompression& GPUrestrict() decompressor = processors.tpcDecompressor; + const ClusterNativeAccess* clusterAccess = decompressor.mClusterNativeAccess; + for (uint32_t i = get_global_id(0); i < GPUCA_NSLICES * GPUCA_ROW_COUNT; i += get_global_size(0)) { + uint32_t slice = i / GPUCA_ROW_COUNT; + uint32_t row = i % GPUCA_ROW_COUNT; + for (uint32_t k = 0; k < clusterAccess->nClusters[slice][row]; k++) { + ClusterNative cl = clusterAccess->clusters[slice][row][k]; + if (isClusterKept(cl, param)) { + decompressor.mNClusterPerSectorRow[i]++; + } + } + } +} + +template <> +GPUdii() void GPUTPCDecompressionUtilKernels::Thread(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread, GPUsharedref() GPUSharedMemory& smem, processorType& processors) +{ + const GPUParam& GPUrestrict() param = processors.param; + GPUTPCDecompression& GPUrestrict() decompressor = processors.tpcDecompressor; + ClusterNative* GPUrestrict() clusterBuffer = decompressor.mNativeClustersBuffer; + const ClusterNativeAccess* clusterAccess = decompressor.mClusterNativeAccess; + const ClusterNativeAccess* outputAccess = processors.ioPtrs.clustersNative; + for (uint32_t i = get_global_id(0); i < GPUCA_NSLICES * GPUCA_ROW_COUNT; i += get_global_size(0)) { + uint32_t slice = i / GPUCA_ROW_COUNT; + uint32_t row = i % GPUCA_ROW_COUNT; + uint32_t count = 0; + for (uint32_t k = 0; k < clusterAccess->nClusters[slice][row]; k++) { + const ClusterNative cl = clusterAccess->clusters[slice][row][k]; + if (isClusterKept(cl, param)) { + clusterBuffer[outputAccess->clusterOffset[slice][row] + count] = cl; + count++; + } + } + } +} + template <> GPUdii() void GPUTPCDecompressionUtilKernels::Thread(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread, GPUsharedref() GPUSharedMemory& smem, processorType& processors) { diff --git a/GPU/GPUTracking/DataCompression/GPUTPCDecompressionKernels.h b/GPU/GPUTracking/DataCompression/GPUTPCDecompressionKernels.h index 622e1fd984fa7..b45af622ebac8 100644 --- a/GPU/GPUTracking/DataCompression/GPUTPCDecompressionKernels.h +++ b/GPU/GPUTracking/DataCompression/GPUTPCDecompressionKernels.h @@ -59,11 +59,15 @@ class GPUTPCDecompressionUtilKernels : public GPUKernelTemplate { public: enum K : int32_t { - sortPerSectorRow = 0, + countFilteredClusters = 0, + storeFilteredClusters = 1, + sortPerSectorRow = 2, }; template GPUd() static void Thread(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread, GPUsharedref() GPUSharedMemory& smem, processorType& GPUrestrict() processors); + + GPUdi() static bool isClusterKept(const o2::tpc::ClusterNative& cl, const GPUParam& GPUrestrict() param); }; } // namespace GPUCA_NAMESPACE::gpu diff --git a/GPU/GPUTracking/DataTypes/GPUNewCalibValues.cxx b/GPU/GPUTracking/DataTypes/GPUNewCalibValues.cxx index f443809d15ef5..f4061fa12873c 100644 --- a/GPU/GPUTracking/DataTypes/GPUNewCalibValues.cxx +++ b/GPU/GPUTracking/DataTypes/GPUNewCalibValues.cxx @@ -19,9 +19,15 @@ using namespace GPUCA_NAMESPACE::gpu; void GPUNewCalibValues::updateFrom(const GPUNewCalibValues* from) { if (from->newSolenoidField) { - solenoidField = from->newSolenoidField; + newSolenoidField = true; + solenoidField = from->solenoidField; } if (from->newContinuousMaxTimeBin) { + newContinuousMaxTimeBin = true; continuousMaxTimeBin = from->continuousMaxTimeBin; } + if (from->newTPCTimeBinCut) { + newTPCTimeBinCut = true; + tpcTimeBinCut = from->tpcTimeBinCut; + } } diff --git a/GPU/GPUTracking/DataTypes/GPUNewCalibValues.h b/GPU/GPUTracking/DataTypes/GPUNewCalibValues.h index 802306e996553..5d5a31785928c 100644 --- a/GPU/GPUTracking/DataTypes/GPUNewCalibValues.h +++ b/GPU/GPUTracking/DataTypes/GPUNewCalibValues.h @@ -25,8 +25,10 @@ namespace gpu struct GPUNewCalibValues { bool newSolenoidField = false; bool newContinuousMaxTimeBin = false; + bool newTPCTimeBinCut = false; float solenoidField = 0.f; uint32_t continuousMaxTimeBin = 0; + int32_t tpcTimeBinCut = 0; void updateFrom(const GPUNewCalibValues* from); }; diff --git a/GPU/GPUTracking/DataTypes/GPUSettings.h b/GPU/GPUTracking/DataTypes/GPUSettings.h index 69bfb15e3f4b0..b967a7ce42620 100644 --- a/GPU/GPUTracking/DataTypes/GPUSettings.h +++ b/GPU/GPUTracking/DataTypes/GPUSettings.h @@ -60,6 +60,7 @@ struct GPUSettingsGRP { int32_t grpContinuousMaxTimeBin = -2; // 0 for triggered events, -1 for automatic setting, -2 invalid default int32_t needsClusterer = 0; // Set to true if the data requires the clusterizer int32_t doCompClusterDecode = 0; // Set to true if the data contains compressed TPC clusters + int32_t tpcCutTimeBin = 0; // Cut TPC clusters and digits >= this cut }; // Parameters of the current time frame diff --git a/GPU/GPUTracking/Definitions/GPUDefGPUParameters.h b/GPU/GPUTracking/Definitions/GPUDefGPUParameters.h index 970e1b2926853..3852d37f6facf 100644 --- a/GPU/GPUTracking/Definitions/GPUDefGPUParameters.h +++ b/GPU/GPUTracking/Definitions/GPUDefGPUParameters.h @@ -344,6 +344,12 @@ #endif #ifndef GPUCA_LB_GPUTPCDecompressionUtilKernels_sortPerSectorRow #define GPUCA_LB_GPUTPCDecompressionUtilKernels_sortPerSectorRow 256 + #endif + #ifndef GPUCA_LB_GPUTPCDecompressionUtilKernels_countFilteredClusters + #define GPUCA_LB_GPUTPCDecompressionUtilKernels_countFilteredClusters 256 + #endif + #ifndef GPUCA_LB_GPUTPCDecompressionUtilKernels_storeFilteredClusters + #define GPUCA_LB_GPUTPCDecompressionUtilKernels_storeFilteredClusters 256 #endif #ifndef GPUCA_LB_GPUTPCCFDecodeZS #define GPUCA_LB_GPUTPCCFDecodeZS 128, 4 diff --git a/GPU/GPUTracking/Definitions/GPUSettingsList.h b/GPU/GPUTracking/Definitions/GPUSettingsList.h index d5494d04930f5..07cd320140909 100644 --- a/GPU/GPUTracking/Definitions/GPUSettingsList.h +++ b/GPU/GPUTracking/Definitions/GPUSettingsList.h @@ -181,6 +181,11 @@ AddOptionRTC(pileupBwdNBC, uint8_t, 80, "", 0, "Pre-trigger Pile-up integration AddHelp("help", 'h') EndConfig() +// Dynamic settings, must NOT use AddOptionRTC(...) !!! +BeginSubConfig(GPUSettingsRecDynamic, dyn, configStandalone.rec, "RECDYN", 0, "Reconstruction settings", rec_dyn) +AddHelp("help", 'h') +EndConfig() + BeginSubConfig(GPUSettingsRec, rec, configStandalone, "REC", 0, "Reconstruction settings", rec) AddOptionRTC(maxTrackQPtB5, float, 1.f / GPUCA_MIN_TRACK_PTB5_DEFAULT, "", 0, "required max Q/Pt (==min Pt) of tracks") AddOptionRTC(nonConsecutiveIDs, int8_t, false, "", 0, "Non-consecutive cluster IDs as in HLT, disables features that need access to slice data in TPC merger") @@ -193,6 +198,7 @@ AddOptionRTC(trackingRefitGPUModel, int8_t, 1, "", 0, "Use GPU track model for t AddCustomCPP(void SetMinTrackPtB5(float v) { maxTrackQPtB5 = v > 0.001f ? (1.f / v) : (1.f / 0.001f); }) AddSubConfig(GPUSettingsRecTPC, tpc) AddSubConfig(GPUSettingsRecTRD, trd) +AddSubConfig(GPUSettingsRecDynamic, dyn) AddHelp("help", 'h') EndConfig() @@ -530,6 +536,7 @@ AddOption(solenoidBzNominalGPU, float, -1e6f, "", 0, "Field strength of solenoid AddOption(constBz, bool, false, "", 0, "force constant Bz for tests") AddOption(setMaxTimeBin, int32_t, -2, "", 0, "maximum time bin of continuous data, 0 for triggered events, -1 for automatic continuous mode, -2 for automatic continuous / triggered") AddOption(overrideNHbfPerTF, int32_t, 0, "", 0, "Overrides the number of HBF per TF if != 0") +AddOption(overrideTPCTimeBinCur, int32_t, 0, "", 0, "Overrides TPC time bin cut if > 0") AddOption(deviceType, std::string, "CPU", "", 0, "Device type, CPU | CUDA | HIP | OCL1 | OCL2") AddOption(forceDeviceType, bool, true, "", 0, "force device type, otherwise allows fall-back to CPU") AddOption(synchronousProcessing, bool, false, "", 0, "Apply performance shortcuts for synchronous processing, disable unneeded steps") diff --git a/GPU/GPUTracking/GPUTrackingLinkDef_O2_DataTypes.h b/GPU/GPUTracking/GPUTrackingLinkDef_O2_DataTypes.h index 16f9b769123f7..6ed4e036c6597 100644 --- a/GPU/GPUTracking/GPUTrackingLinkDef_O2_DataTypes.h +++ b/GPU/GPUTracking/GPUTrackingLinkDef_O2_DataTypes.h @@ -28,6 +28,7 @@ #pragma link C++ class o2::gpu::GPUConfigurableParamGPUSettingsRec + ; #pragma link C++ class o2::gpu::GPUConfigurableParamGPUSettingsRecTPC + ; #pragma link C++ class o2::gpu::GPUConfigurableParamGPUSettingsRecTRD + ; +#pragma link C++ class o2::gpu::GPUConfigurableParamGPUSettingsRecDynamic + ; #pragma link C++ class o2::gpu::GPUConfigurableParamGPUSettingsProcessing + ; #pragma link C++ class o2::gpu::GPUConfigurableParamGPUSettingsProcessingParam + ; #pragma link C++ class o2::gpu::GPUConfigurableParamGPUSettingsProcessingRTC + ; diff --git a/GPU/GPUTracking/Global/GPUChainTracking.cxx b/GPU/GPUTracking/Global/GPUChainTracking.cxx index 8c2599604387b..ff476716febe8 100644 --- a/GPU/GPUTracking/Global/GPUChainTracking.cxx +++ b/GPU/GPUTracking/Global/GPUChainTracking.cxx @@ -633,7 +633,7 @@ int32_t GPUChainTracking::DoQueuedUpdates(int32_t stream, bool updateSlave) const GPUSettingsProcessing* p = nullptr; std::lock_guard lk(mMutexUpdateCalib); if (mUpdateNewCalibObjects) { - if (mNewCalibValues->newSolenoidField || mNewCalibValues->newContinuousMaxTimeBin) { + if (mNewCalibValues->newSolenoidField || mNewCalibValues->newContinuousMaxTimeBin || mNewCalibValues->newTPCTimeBinCut) { grp = std::make_unique(mRec->GetGRPSettings()); if (mNewCalibValues->newSolenoidField) { grp->solenoidBzNominalGPU = mNewCalibValues->solenoidField; @@ -641,6 +641,9 @@ int32_t GPUChainTracking::DoQueuedUpdates(int32_t stream, bool updateSlave) if (mNewCalibValues->newContinuousMaxTimeBin) { grp->grpContinuousMaxTimeBin = mNewCalibValues->continuousMaxTimeBin; } + if (mNewCalibValues->newTPCTimeBinCut) { + grp->tpcCutTimeBin = mNewCalibValues->tpcTimeBinCut; + } } } if (GetProcessingSettings().tpcDownscaledEdx != 0) { @@ -921,33 +924,6 @@ int32_t GPUChainTracking::FinalizePipelinedProcessing() return RunChainFinalize(); } -int32_t GPUChainTracking::HelperReadEvent(int32_t iSlice, int32_t threadId, GPUReconstructionHelpers::helperParam* par) { return ReadEvent(iSlice, threadId); } - -int32_t GPUChainTracking::HelperOutput(int32_t iSlice, int32_t threadId, GPUReconstructionHelpers::helperParam* par) -{ - if (param().rec.tpc.globalTracking) { - uint32_t tmpSlice = GPUTPCGlobalTracking::GlobalTrackingSliceOrder(iSlice); - uint32_t sliceLeft, sliceRight; - GPUTPCGlobalTracking::GlobalTrackingSliceLeftRight(tmpSlice, sliceLeft, sliceRight); - - while (mSliceSelectorReady < (int32_t)tmpSlice || mSliceSelectorReady < (int32_t)sliceLeft || mSliceSelectorReady < (int32_t)sliceRight) { - if (par->reset) { - return 1; - } - } - GlobalTracking(tmpSlice, 0); - WriteOutput(tmpSlice, 0); - } else { - while (mSliceSelectorReady < iSlice) { - if (par->reset) { - return 1; - } - } - WriteOutput(iSlice, threadId); - } - return 0; -} - int32_t GPUChainTracking::CheckErrorCodes(bool cpuOnly, bool forceShowErrors, std::vector>* fillErrors) { int32_t retVal = 0; diff --git a/GPU/GPUTracking/Global/GPUChainTrackingClusterizer.cxx b/GPU/GPUTracking/Global/GPUChainTrackingClusterizer.cxx index ae240181eba65..4bc0ee4e91ff1 100644 --- a/GPU/GPUTracking/Global/GPUChainTrackingClusterizer.cxx +++ b/GPU/GPUTracking/Global/GPUChainTrackingClusterizer.cxx @@ -18,6 +18,7 @@ #include "GPUO2DataTypes.h" #include "GPUMemorySizeScalers.h" #include "GPUTrackingInputProvider.h" +#include "GPUNewCalibValues.h" #include #ifdef GPUCA_O2_LIB @@ -457,6 +458,7 @@ std::pair GPUChainTracking::RunTPCClusterizer_transferZS(int int32_t GPUChainTracking::RunTPCClusterizer_prepare(bool restorePointers) { + bool doGPU = mRec->GetRecoStepsGPU() & GPUDataTypes::RecoStep::TPCClusterFinding; if (restorePointers) { for (uint32_t iSlice = 0; iSlice < NSLICES; iSlice++) { processors()->tpcClusterer[iSlice].mPzsOffsets = mCFContext->ptrSave[iSlice].zsOffsetHost; @@ -512,7 +514,7 @@ int32_t GPUChainTracking::RunTPCClusterizer_prepare(bool restorePointers) uint32_t threshold = 40000000; uint32_t nDigitsScaled = nDigitsBase > threshold ? nDigitsBase : std::min((threshold + nDigitsBase) / 2, 2 * nDigitsBase); processors()->tpcClusterer[iSlice].SetNMaxDigits(processors()->tpcClusterer[iSlice].mPmemory->counters.nDigits, mCFContext->nPagesFragmentMax, nDigitsScaled, mCFContext->nDigitsEndpointMax[iSlice]); - if (mRec->IsGPU()) { + if (doGPU) { processorsShadow()->tpcClusterer[iSlice].SetNMaxDigits(processors()->tpcClusterer[iSlice].mPmemory->counters.nDigits, mCFContext->nPagesFragmentMax, nDigitsScaled, mCFContext->nDigitsEndpointMax[iSlice]); } if (mPipelineNotifyCtx && GetProcessingSettings().doublePipelineClusterizer) { @@ -572,13 +574,14 @@ int32_t GPUChainTracking::RunTPCClusterizer(bool synchronizeOutput) return ForwardTPCDigits(); } #ifdef GPUCA_TPC_GEOMETRY_O2 + int32_t tpcTimeBinCut = mUpdateNewCalibObjects && mNewCalibValues->newTPCTimeBinCut ? mNewCalibValues->tpcTimeBinCut : param().tpcCutTimeBin; mRec->PushNonPersistentMemory(qStr2Tag("TPCCLUST")); const auto& threadContext = GetThreadContext(); const bool doGPU = GetRecoStepsGPU() & RecoStep::TPCClusterFinding; if (RunTPCClusterizer_prepare(mPipelineNotifyCtx && GetProcessingSettings().doublePipelineClusterizer)) { return 1; } - if (GetProcessingSettings().ompAutoNThreads && !mRec->IsGPU()) { + if (GetProcessingSettings().ompAutoNThreads && !doGPU) { mRec->SetNOMPThreads(mRec->MemoryScalers()->nTPCdigits / 20000); } @@ -764,6 +767,7 @@ int32_t GPUChainTracking::RunTPCClusterizer(bool synchronizeOutput) : 0; uint32_t nBlocks = doGPU ? clusterer.mPmemory->counters.nPagesSubslice : GPUTrackingInOutZS::NENDPOINTS; + (void)tpcTimeBinCut; // TODO: To be used in decoding kernels switch (mCFContext->zsVersion) { default: GPUFatal("Data with invalid TPC ZS mode (%d) received", mCFContext->zsVersion); diff --git a/GPU/GPUTracking/Global/GPUChainTrackingCompression.cxx b/GPU/GPUTracking/Global/GPUChainTrackingCompression.cxx index 98109447de034..01e4d011d08b9 100644 --- a/GPU/GPUTracking/Global/GPUChainTrackingCompression.cxx +++ b/GPU/GPUTracking/Global/GPUChainTrackingCompression.cxx @@ -219,14 +219,15 @@ int32_t GPUChainTracking::RunTPCDecompression() return ((tmpBuffer = std::make_unique(size))).get(); }; auto& decompressTimer = getTimer("TPCDecompression", 0); - auto allocatorUse = GetProcessingSettings().tpcApplyCFCutsAtDecoding ? std::function{allocatorTmp} : std::function{allocatorFinal}; + bool runFiltering = GetProcessingSettings().tpcApplyCFCutsAtDecoding; + auto allocatorUse = runFiltering ? std::function{allocatorTmp} : std::function{allocatorFinal}; decompressTimer.Start(); if (decomp.decompress(mIOPtrs.tpcCompressedClusters, *mClusterNativeAccess, allocatorUse, param(), GetProcessingSettings().deterministicGPUReconstruction)) { GPUError("Error decompressing clusters"); return 1; } - if (GetProcessingSettings().tpcApplyCFCutsAtDecoding) { - RunTPCClusterFilter(mClusterNativeAccess.get(), allocatorFinal, true); + if (runFiltering) { + RunTPCClusterFilter(mClusterNativeAccess.get(), allocatorFinal, GetProcessingSettings().tpcApplyCFCutsAtDecoding); } decompressTimer.Stop(); mIOPtrs.clustersNative = mClusterNativeAccess.get(); @@ -245,6 +246,7 @@ int32_t GPUChainTracking::RunTPCDecompression() mRec->PushNonPersistentMemory(qStr2Tag("TPCDCMPR")); RecoStep myStep = RecoStep::TPCDecompression; bool doGPU = GetRecoStepsGPU() & RecoStep::TPCDecompression; + bool runFiltering = param().tpcCutTimeBin > 0; GPUTPCDecompression& Decompressor = processors()->tpcDecompressor; GPUTPCDecompression& DecompressorShadow = doGPU ? processorsShadow()->tpcDecompressor : Decompressor; const auto& threadContext = GetThreadContext(); @@ -252,42 +254,19 @@ int32_t GPUChainTracking::RunTPCDecompression() CompressedClusters& inputGPU = Decompressor.mInputGPU; CompressedClusters& inputGPUShadow = DecompressorShadow.mInputGPU; + if (cmprClsHost.nTracks && cmprClsHost.solenoidBz != -1e6f && cmprClsHost.solenoidBz != param().bzkG) { + throw std::runtime_error("Configured solenoid Bz does not match value used for track model encoding"); + } + if (cmprClsHost.nTracks && cmprClsHost.maxTimeBin != -1e6 && cmprClsHost.maxTimeBin != param().continuousMaxTimeBin) { + throw std::runtime_error("Configured max time bin does not match value used for track model encoding"); + } + int32_t inputStream = 0; int32_t unattachedStream = mRec->NStreams() - 1; - inputGPU.nAttachedClusters = cmprClsHost.nAttachedClusters; - inputGPU.nUnattachedClusters = cmprClsHost.nUnattachedClusters; - inputGPU.nTracks = cmprClsHost.nTracks; - inputGPU.nAttachedClustersReduced = inputGPU.nAttachedClusters - inputGPU.nTracks; - inputGPU.nSliceRows = NSLICES * GPUCA_ROW_COUNT; - inputGPU.nComppressionModes = param().rec.tpc.compressionTypeMask; - inputGPU.solenoidBz = param().bzkG; - inputGPU.maxTimeBin = param().continuousMaxTimeBin; + inputGPU = cmprClsHost; SetupGPUProcessor(&Decompressor, true); WriteToConstantMemory(myStep, (char*)&processors()->tpcDecompressor - (char*)processors(), &DecompressorShadow, sizeof(DecompressorShadow), inputStream); - - inputGPU.nTrackClusters = cmprClsHost.nTrackClusters; - inputGPU.qTotU = cmprClsHost.qTotU; - inputGPU.qMaxU = cmprClsHost.qMaxU; - inputGPU.flagsU = cmprClsHost.flagsU; - inputGPU.padDiffU = cmprClsHost.padDiffU; - inputGPU.timeDiffU = cmprClsHost.timeDiffU; - inputGPU.sigmaPadU = cmprClsHost.sigmaPadU; - inputGPU.sigmaTimeU = cmprClsHost.sigmaTimeU; - inputGPU.nSliceRowClusters = cmprClsHost.nSliceRowClusters; - inputGPU.qTotA = cmprClsHost.qTotA; - inputGPU.qMaxA = cmprClsHost.qMaxA; - inputGPU.flagsA = cmprClsHost.flagsA; - inputGPU.rowDiffA = cmprClsHost.rowDiffA; - inputGPU.sliceLegDiffA = cmprClsHost.sliceLegDiffA; - inputGPU.padResA = cmprClsHost.padResA; - inputGPU.timeResA = cmprClsHost.timeResA; - inputGPU.sigmaPadA = cmprClsHost.sigmaPadA; - inputGPU.sigmaTimeA = cmprClsHost.sigmaTimeA; - inputGPU.qPtA = cmprClsHost.qPtA; - inputGPU.rowA = cmprClsHost.rowA; - inputGPU.sliceA = cmprClsHost.sliceA; - inputGPU.timeA = cmprClsHost.timeA; - inputGPU.padA = cmprClsHost.padA; + inputGPU = cmprClsHost; bool toGPU = true; runKernel({GetGridAutoStep(inputStream, RecoStep::TPCDecompression), krnlRunRangeNone, &mEvents->init}, DecompressorShadow.mNativeClustersIndex, NSLICES * GPUCA_ROW_COUNT * sizeof(DecompressorShadow.mNativeClustersIndex[0])); @@ -329,12 +308,6 @@ int32_t GPUChainTracking::RunTPCDecompression() GPUMemCpy(myStep, inputGPUShadow.sigmaPadU, cmprClsHost.sigmaPadU, cmprClsHost.nUnattachedClusters * sizeof(cmprClsHost.sigmaPadU[0]), unattachedStream, toGPU); GPUMemCpy(myStep, inputGPUShadow.sigmaTimeU, cmprClsHost.sigmaTimeU, cmprClsHost.nUnattachedClusters * sizeof(cmprClsHost.sigmaTimeU[0]), unattachedStream, toGPU); - mInputsHost->mNClusterNative = mInputsShadow->mNClusterNative = cmprClsHost.nAttachedClusters + cmprClsHost.nUnattachedClusters; - AllocateRegisteredMemory(mInputsHost->mResourceClusterNativeOutput, mSubOutputControls[GPUTrackingOutputs::getIndex(&GPUTrackingOutputs::clustersNative)]); - AllocateRegisteredMemory(mInputsHost->mResourceClusterNativeBuffer); - DecompressorShadow.mNativeClustersBuffer = mInputsShadow->mPclusterNativeBuffer; - Decompressor.mNativeClustersBuffer = mInputsHost->mPclusterNativeOutput; - WriteToConstantMemory(myStep, (char*)&processors()->tpcDecompressor - (char*)processors(), &DecompressorShadow, sizeof(DecompressorShadow), inputStream); TransferMemoryResourceLinkToHost(RecoStep::TPCDecompression, Decompressor.mResourceTmpIndexes, inputStream, nullptr, mEvents->stream, nStreams); SynchronizeStream(inputStream); uint32_t offset = 0; @@ -353,27 +326,83 @@ int32_t GPUChainTracking::RunTPCDecompression() if (decodedAttachedClusters != cmprClsHost.nAttachedClusters) { GPUWarning("%u / %u clusters failed track model decoding (%f %%)", cmprClsHost.nAttachedClusters - decodedAttachedClusters, cmprClsHost.nAttachedClusters, 100.f * (float)(cmprClsHost.nAttachedClusters - decodedAttachedClusters) / (float)cmprClsHost.nAttachedClusters); } - if (doGPU) { - mClusterNativeAccess->clustersLinear = mInputsShadow->mPclusterNativeBuffer; + if (runFiltering) { // If filtering, allocate a temporary buffer and cluster native access in decompressor context + Decompressor.mNClusterNativeBeforeFiltering = DecompressorShadow.mNClusterNativeBeforeFiltering = decodedAttachedClusters + cmprClsHost.nUnattachedClusters; + AllocateRegisteredMemory(Decompressor.mResourceTmpBufferBeforeFiltering); + AllocateRegisteredMemory(Decompressor.mResourceClusterNativeAccess); + mClusterNativeAccess->clustersLinear = DecompressorShadow.mNativeClustersBuffer; + mClusterNativeAccess->setOffsetPtrs(); + *Decompressor.mClusterNativeAccess = *mClusterNativeAccess; + WriteToConstantMemory(myStep, (char*)&processors()->tpcDecompressor - (char*)processors(), &DecompressorShadow, sizeof(DecompressorShadow), inputStream); + TransferMemoryResourceLinkToGPU(RecoStep::TPCDecompression, Decompressor.mResourceClusterNativeAccess, inputStream, &mEvents->single); + } else { // If not filtering, directly allocate the final buffers + mInputsHost->mNClusterNative = mInputsShadow->mNClusterNative = cmprClsHost.nAttachedClusters + cmprClsHost.nUnattachedClusters; + AllocateRegisteredMemory(mInputsHost->mResourceClusterNativeOutput, mSubOutputControls[GPUTrackingOutputs::getIndex(&GPUTrackingOutputs::clustersNative)]); + AllocateRegisteredMemory(mInputsHost->mResourceClusterNativeBuffer); + DecompressorShadow.mNativeClustersBuffer = mInputsShadow->mPclusterNativeBuffer; + Decompressor.mNativeClustersBuffer = mInputsHost->mPclusterNativeOutput; + DecompressorShadow.mClusterNativeAccess = mInputsShadow->mPclusterNativeAccess; + Decompressor.mClusterNativeAccess = mInputsHost->mPclusterNativeAccess; + WriteToConstantMemory(myStep, (char*)&processors()->tpcDecompressor - (char*)processors(), &DecompressorShadow, sizeof(DecompressorShadow), inputStream); + if (doGPU) { + mClusterNativeAccess->clustersLinear = mInputsShadow->mPclusterNativeBuffer; + mClusterNativeAccess->setOffsetPtrs(); + *mInputsHost->mPclusterNativeAccess = *mClusterNativeAccess; + processorsShadow()->ioPtrs.clustersNative = mInputsShadow->mPclusterNativeAccess; + WriteToConstantMemory(RecoStep::TPCDecompression, (char*)&processors()->ioPtrs - (char*)processors(), &processorsShadow()->ioPtrs, sizeof(processorsShadow()->ioPtrs), inputStream); + TransferMemoryResourceLinkToGPU(RecoStep::TPCDecompression, mInputsHost->mResourceClusterNativeAccess, inputStream, &mEvents->single); + } + mIOPtrs.clustersNative = mClusterNativeAccess.get(); + mClusterNativeAccess->clustersLinear = mInputsHost->mPclusterNativeOutput; mClusterNativeAccess->setOffsetPtrs(); *mInputsHost->mPclusterNativeAccess = *mClusterNativeAccess; - processorsShadow()->ioPtrs.clustersNative = mInputsShadow->mPclusterNativeAccess; - WriteToConstantMemory(RecoStep::TPCDecompression, (char*)&processors()->ioPtrs - (char*)processors(), &processorsShadow()->ioPtrs, sizeof(processorsShadow()->ioPtrs), inputStream); - TransferMemoryResourceLinkToGPU(RecoStep::TPCDecompression, mInputsHost->mResourceClusterNativeAccess, inputStream, &mEvents->single); } - mIOPtrs.clustersNative = mClusterNativeAccess.get(); - mClusterNativeAccess->clustersLinear = mInputsHost->mPclusterNativeOutput; - mClusterNativeAccess->setOffsetPtrs(); uint32_t batchSize = doGPU ? 6 : NSLICES; for (uint32_t iSlice = 0; iSlice < NSLICES; iSlice = iSlice + batchSize) { int32_t iStream = (iSlice / batchSize) % mRec->NStreams(); runKernel({GetGridAuto(iStream), krnlRunRangeNone, {nullptr, &mEvents->single}}, iSlice, batchSize); uint32_t copySize = std::accumulate(mClusterNativeAccess->nClustersSector + iSlice, mClusterNativeAccess->nClustersSector + iSlice + batchSize, 0u); - GPUMemCpy(RecoStep::TPCDecompression, mInputsHost->mPclusterNativeOutput + mClusterNativeAccess->clusterOffset[iSlice][0], DecompressorShadow.mNativeClustersBuffer + mClusterNativeAccess->clusterOffset[iSlice][0], sizeof(Decompressor.mNativeClustersBuffer[0]) * copySize, iStream, false); + if (!runFiltering) { + GPUMemCpy(RecoStep::TPCDecompression, mInputsHost->mPclusterNativeOutput + mClusterNativeAccess->clusterOffset[iSlice][0], DecompressorShadow.mNativeClustersBuffer + mClusterNativeAccess->clusterOffset[iSlice][0], sizeof(Decompressor.mNativeClustersBuffer[0]) * copySize, iStream, false); + } } SynchronizeGPU(); + if (runFiltering) { // If filtering is applied, count how many clusters will remain after filtering and allocate final buffers accordingly + AllocateRegisteredMemory(Decompressor.mResourceNClusterPerSectorRow); + WriteToConstantMemory(myStep, (char*)&processors()->tpcDecompressor - (char*)processors(), &DecompressorShadow, sizeof(DecompressorShadow), unattachedStream); + runKernel({GetGridAutoStep(unattachedStream, RecoStep::TPCDecompression), krnlRunRangeNone}, DecompressorShadow.mNClusterPerSectorRow, NSLICES * GPUCA_ROW_COUNT * sizeof(DecompressorShadow.mNClusterPerSectorRow[0])); + runKernel(GetGridAutoStep(unattachedStream, RecoStep::TPCDecompression)); + TransferMemoryResourceLinkToHost(RecoStep::TPCDecompression, Decompressor.mResourceNClusterPerSectorRow, unattachedStream); + SynchronizeStream(unattachedStream); + uint32_t nClustersFinal = std::accumulate(Decompressor.mNClusterPerSectorRow, Decompressor.mNClusterPerSectorRow + inputGPU.nSliceRows, 0u); + mInputsHost->mNClusterNative = mInputsShadow->mNClusterNative = nClustersFinal; + AllocateRegisteredMemory(mInputsHost->mResourceClusterNativeOutput, mSubOutputControls[GPUTrackingOutputs::getIndex(&GPUTrackingOutputs::clustersNative)]); + AllocateRegisteredMemory(mInputsHost->mResourceClusterNativeBuffer); + DecompressorShadow.mNativeClustersBuffer = mInputsShadow->mPclusterNativeBuffer; + Decompressor.mNativeClustersBuffer = mInputsHost->mPclusterNativeOutput; + WriteToConstantMemory(myStep, (char*)&processors()->tpcDecompressor - (char*)processors(), &DecompressorShadow, sizeof(DecompressorShadow), unattachedStream); + for (uint32_t i = 0; i < NSLICES; i++) { + for (uint32_t j = 0; j < GPUCA_ROW_COUNT; j++) { + mClusterNativeAccess->nClusters[i][j] = Decompressor.mNClusterPerSectorRow[i * GPUCA_ROW_COUNT + j]; + } + } + if (doGPU) { + mClusterNativeAccess->clustersLinear = mInputsShadow->mPclusterNativeBuffer; + mClusterNativeAccess->setOffsetPtrs(); + *mInputsHost->mPclusterNativeAccess = *mClusterNativeAccess; + processorsShadow()->ioPtrs.clustersNative = mInputsShadow->mPclusterNativeAccess; + WriteToConstantMemory(RecoStep::TPCDecompression, (char*)&processors()->ioPtrs - (char*)processors(), &processorsShadow()->ioPtrs, sizeof(processorsShadow()->ioPtrs), unattachedStream); + TransferMemoryResourceLinkToGPU(RecoStep::TPCDecompression, mInputsHost->mResourceClusterNativeAccess, unattachedStream); + } + mIOPtrs.clustersNative = mClusterNativeAccess.get(); + mClusterNativeAccess->clustersLinear = mInputsHost->mPclusterNativeOutput; + mClusterNativeAccess->setOffsetPtrs(); + runKernel(GetGridAutoStep(unattachedStream, RecoStep::TPCDecompression)); + GPUMemCpy(RecoStep::TPCDecompression, mInputsHost->mPclusterNativeOutput, DecompressorShadow.mNativeClustersBuffer, sizeof(Decompressor.mNativeClustersBuffer[0]) * nClustersFinal, unattachedStream, false); + SynchronizeStream(unattachedStream); + } if (GetProcessingSettings().deterministicGPUReconstruction || GetProcessingSettings().debugLevel >= 4) { runKernel(GetGridAutoStep(unattachedStream, RecoStep::TPCDecompression)); const ClusterNativeAccess* decoded = mIOPtrs.clustersNative; @@ -386,6 +415,7 @@ int32_t GPUChainTracking::RunTPCDecompression() } } } + SynchronizeStream(unattachedStream); } mRec->PopNonPersistentMemory(RecoStep::TPCDecompression, qStr2Tag("TPCDCMPR")); } diff --git a/GPU/GPUTracking/Global/GPUChainTrackingDebugAndProfiling.cxx b/GPU/GPUTracking/Global/GPUChainTrackingDebugAndProfiling.cxx index f8a64e9d4faaa..7d4a3420995ad 100644 --- a/GPU/GPUTracking/Global/GPUChainTrackingDebugAndProfiling.cxx +++ b/GPU/GPUTracking/Global/GPUChainTrackingDebugAndProfiling.cxx @@ -315,6 +315,9 @@ void GPUChainTracking::RunTPCClusterFilter(o2::tpc::ClusterNativeAccess* cluster keep = keep && cl.qTot > param().rec.tpc.cfQTotCutoff && cl.qMax > param().rec.tpc.cfQMaxCutoff; keep = keep && (!(cl.getFlags() & o2::tpc::ClusterNative::flagSingle) || ((cl.sigmaPadPacked || cl.qMax > param().rec.tpc.cfQMaxCutoffSinglePad) && (cl.sigmaTimePacked || cl.qMax > param().rec.tpc.cfQMaxCutoffSingleTime))); } + if (param().tpcCutTimeBin > 0) { + keep = keep && cl.getTime() < param().tpcCutTimeBin; + } keep = keep && (!GetProcessingSettings().tpcApplyDebugClusterFilter || clusterFilter.filter(iSector, iRow, cl)); if (iPhase && keep) { outputBuffer[countTotal] = cl; diff --git a/GPU/GPUTracking/Global/GPUChainTrackingSliceTracker.cxx b/GPU/GPUTracking/Global/GPUChainTrackingSliceTracker.cxx index 8db15fb1aef7e..62c93bcb1bfb5 100644 --- a/GPU/GPUTracking/Global/GPUChainTrackingSliceTracker.cxx +++ b/GPU/GPUTracking/Global/GPUChainTrackingSliceTracker.cxx @@ -532,3 +532,30 @@ void GPUChainTracking::WriteOutput(int32_t iSlice, int32_t threadId) GPUInfo("Finished WriteOutput for slice %d on thread %d\n", iSlice, threadId); } } + +int32_t GPUChainTracking::HelperReadEvent(int32_t iSlice, int32_t threadId, GPUReconstructionHelpers::helperParam* par) { return ReadEvent(iSlice, threadId); } + +int32_t GPUChainTracking::HelperOutput(int32_t iSlice, int32_t threadId, GPUReconstructionHelpers::helperParam* par) +{ + if (param().rec.tpc.globalTracking) { + uint32_t tmpSlice = GPUTPCGlobalTracking::GlobalTrackingSliceOrder(iSlice); + uint32_t sliceLeft, sliceRight; + GPUTPCGlobalTracking::GlobalTrackingSliceLeftRight(tmpSlice, sliceLeft, sliceRight); + + while (mSliceSelectorReady < (int32_t)tmpSlice || mSliceSelectorReady < (int32_t)sliceLeft || mSliceSelectorReady < (int32_t)sliceRight) { + if (par->reset) { + return 1; + } + } + GlobalTracking(tmpSlice, 0); + WriteOutput(tmpSlice, 0); + } else { + while (mSliceSelectorReady < iSlice) { + if (par->reset) { + return 1; + } + } + WriteOutput(iSlice, threadId); + } + return 0; +} diff --git a/GPU/GPUTracking/Merger/GPUTPCGMTrackParam.cxx b/GPU/GPUTracking/Merger/GPUTPCGMTrackParam.cxx index d817278404534..74cc12e9bbd9a 100644 --- a/GPU/GPUTracking/Merger/GPUTPCGMTrackParam.cxx +++ b/GPU/GPUTracking/Merger/GPUTPCGMTrackParam.cxx @@ -370,7 +370,7 @@ GPUd() bool GPUTPCGMTrackParam::Fit(GPUTPCGMMerger* GPUrestrict() merger, int32_ CADEBUG(printf("Reinit linearization\n")); prop.SetTrack(this, prop.GetAlpha()); } - if (param.par.dodEdx && param.dodEdxDownscaled && iWay == nWays - 1 && cluster.leg == clusters[maxN - 1].leg && !(clusterState & GPUTPCGMMergedTrackHit::flagEdge)) { + if (param.par.dodEdx && param.dodEdxDownscaled && iWay == nWays - 1 && cluster.leg == clusters[maxN - 1].leg && !(clusterState & GPUTPCGMMergedTrackHit::flagEdge)) { // TODO: Costimize flag to remove, and option to remove double-clusters float qtot = 0, qmax = 0, pad = 0, relTime = 0; const int32_t clusterCount = (ihit - ihitMergeFirst) * wayDirection + 1; for (int32_t iTmp = ihitMergeFirst; iTmp != ihit + wayDirection; iTmp += wayDirection) { @@ -384,7 +384,7 @@ GPUd() bool GPUTPCGMTrackParam::Fit(GPUTPCGMMerger* GPUrestrict() merger, int32_ relTime += cl.getTime(); } } - qtot /= clusterCount; + qtot /= clusterCount; // TODO: Weighted Average pad /= clusterCount; relTime /= clusterCount; relTime = relTime - CAMath::Round(relTime); @@ -528,7 +528,7 @@ GPUd() int32_t GPUTPCGMTrackParam::MergeDoubleRowClusters(int32_t& ihit, int32_t } } else { CADEBUG(printf("\t\tMerging hit row %d X %f Y %f Z %f (dy %f, dz %f, chiY %f, chiZ %f)\n", clusters[ihit].row, clx, cly, clz, dy, dz, sqrtf(maxDistY), sqrtf(maxDistZ))); - xx += clx * clamp; + xx += clx * clamp; // TODO: Weight in pad/time instead of XYZ yy += cly * clamp; zz += clz * clamp; clusterState |= clusters[ihit].state; diff --git a/GPU/GPUTracking/SliceTracker/GPUTPCSliceData.cxx b/GPU/GPUTracking/SliceTracker/GPUTPCSliceData.cxx index 6908bc326a535..6c456a28918ab 100644 --- a/GPU/GPUTracking/SliceTracker/GPUTPCSliceData.cxx +++ b/GPU/GPUTracking/SliceTracker/GPUTPCSliceData.cxx @@ -265,22 +265,20 @@ GPUdii() int32_t GPUTPCSliceData::InitFromClusterData(int32_t nBlocks, int32_t n for (uint32_t i = iThread; i < NumberOfClusters; i += nThreads) { UpdateMinMaxYZ(yMin, yMax, zMin, zMax, YZData[RowOffset + i].x, YZData[RowOffset + i].y); } + } else if (mem->param.par.earlyTpcTransform) { // Early transform case with ClusterNative present + for (uint32_t i = iThread; i < NumberOfClusters; i += nThreads) { + float2 tmp; + tmp.x = mClusterData[RowOffset + i].y; + tmp.y = mClusterData[RowOffset + i].z; + UpdateMinMaxYZ(yMin, yMax, zMin, zMax, tmp.x, tmp.y); + YZData[RowOffset + i] = tmp; + } } else { - if (mem->param.par.earlyTpcTransform) { // Early transform case with ClusterNative present - for (uint32_t i = iThread; i < NumberOfClusters; i += nThreads) { - float2 tmp; - tmp.x = mClusterData[RowOffset + i].y; - tmp.y = mClusterData[RowOffset + i].z; - UpdateMinMaxYZ(yMin, yMax, zMin, zMax, tmp.x, tmp.y); - YZData[RowOffset + i] = tmp; - } - } else { - for (uint32_t i = iThread; i < NumberOfClusters; i += nThreads) { - float x, y, z; - GPUTPCConvertImpl::convert(*mem, iSlice, rowIndex, mem->ioPtrs.clustersNative->clusters[iSlice][rowIndex][i].getPad(), mem->ioPtrs.clustersNative->clusters[iSlice][rowIndex][i].getTime(), x, y, z); - UpdateMinMaxYZ(yMin, yMax, zMin, zMax, y, z); - YZData[RowOffset + i] = CAMath::MakeFloat2(y, z); - } + for (uint32_t i = iThread; i < NumberOfClusters; i += nThreads) { + float x, y, z; + GPUTPCConvertImpl::convert(*mem, iSlice, rowIndex, mem->ioPtrs.clustersNative->clusters[iSlice][rowIndex][i].getPad(), mem->ioPtrs.clustersNative->clusters[iSlice][rowIndex][i].getTime(), x, y, z); + UpdateMinMaxYZ(yMin, yMax, zMin, zMax, y, z); + YZData[RowOffset + i] = CAMath::MakeFloat2(y, z); } } diff --git a/GPU/GPUTracking/SliceTracker/GPUTPCTracker.cxx b/GPU/GPUTracking/SliceTracker/GPUTPCTracker.cxx index 7428a4ccbd0ed..84bdc52ab6f46 100644 --- a/GPU/GPUTracking/SliceTracker/GPUTPCTracker.cxx +++ b/GPU/GPUTracking/SliceTracker/GPUTPCTracker.cxx @@ -81,7 +81,7 @@ void* GPUTPCTracker::SetPointersScratch(void* mem) if (mRec->GetProcessingSettings().memoryAllocationStrategy != GPUMemoryResource::ALLOCATION_INDIVIDUAL) { mem = SetPointersTracklets(mem); } - if (mRec->IsGPU()) { + if (mRec->GetRecoStepsGPU() & GPUDataTypes::RecoStep::TPCSliceTracking) { computePointerWithAlignment(mem, mTrackletTmpStartHits, GPUCA_ROW_COUNT * mNMaxRowStartHits); computePointerWithAlignment(mem, mRowStartHitCountOffset, GPUCA_ROW_COUNT); } @@ -164,7 +164,7 @@ void GPUTPCTracker::SetMaxData(const GPUTrackingInOutPointers& io) mNMaxTracks = mRec->MemoryScalers()->NTPCSectorTracks(mData.NumberOfHits()); mNMaxTrackHits = mRec->MemoryScalers()->NTPCSectorTrackHits(mData.NumberOfHits(), mRec->GetProcessingSettings().tpcInputWithClusterRejection); #ifdef GPUCA_SORT_STARTHITS_GPU - if (mRec->IsGPU()) { + if (mRec->GetRecoStepsGPU() & GPUDataTypes::RecoStep::TPCSliceTracking) { if (mNMaxStartHits > mNMaxRowStartHits * GPUCA_ROW_COUNT) { mNMaxStartHits = mNMaxRowStartHits * GPUCA_ROW_COUNT; } diff --git a/GPU/GPUTracking/TPCClusterFinder/ClusterAccumulator.cxx b/GPU/GPUTracking/TPCClusterFinder/ClusterAccumulator.cxx index 8988126f7a15e..e8176ecb60d78 100644 --- a/GPU/GPUTracking/TPCClusterFinder/ClusterAccumulator.cxx +++ b/GPU/GPUTracking/TPCClusterFinder/ClusterAccumulator.cxx @@ -27,6 +27,10 @@ GPUd() bool ClusterAccumulator::toNative(const ChargePos& pos, Charge q, tpc::Cl if (cn.qTot <= param.rec.tpc.cfQTotCutoff) { return false; } + cn.qMax = q; + if (cn.qMax <= param.rec.tpc.cfQMaxCutoff) { + return false; + } if (mTimeMean < param.rec.tpc.clustersShiftTimebinsClusterizer) { return false; } @@ -48,7 +52,6 @@ GPUd() bool ClusterAccumulator::toNative(const ChargePos& pos, Charge q, tpc::Cl flags |= (wasSplitInPad) ? tpc::ClusterNative::flagSplitPad : 0; flags |= (isSingleCluster) ? tpc::ClusterNative::flagSingle : 0; - cn.qMax = q; cn.setTimeFlags(mTimeMean - param.rec.tpc.clustersShiftTimebinsClusterizer, flags); cn.setPad(mPadMean); cn.setSigmaTime(mTimeSigma); diff --git a/GPU/GPUTracking/kernels.cmake b/GPU/GPUTracking/kernels.cmake index b0aed5aba1166..f028c6990f267 100644 --- a/GPU/GPUTracking/kernels.cmake +++ b/GPU/GPUTracking/kernels.cmake @@ -108,6 +108,8 @@ o2_gpu_add_kernel("GPUTPCCompressionGatherKernels, multiBlock" "GPUTPCCom o2_gpu_add_kernel("GPUTPCDecompressionKernels, step0attached" "= TPCDECOMPRESSION" LB simple int32_t trackStart int32_t trackEnd) o2_gpu_add_kernel("GPUTPCDecompressionKernels, step1unattached" "= TPCDECOMPRESSION" LB simple int32_t sliceStart int32_t nSlices) o2_gpu_add_kernel("GPUTPCDecompressionUtilKernels, sortPerSectorRow" "GPUTPCDecompressionKernels" LB simple) +o2_gpu_add_kernel("GPUTPCDecompressionUtilKernels, countFilteredClusters" "GPUTPCDecompressionKernels" LB simple) +o2_gpu_add_kernel("GPUTPCDecompressionUtilKernels, storeFilteredClusters" "GPUTPCDecompressionKernels" LB simple) o2_gpu_add_kernel("GPUTPCCFCheckPadBaseline" "= TPCCLUSTERFINDER" LB single) o2_gpu_add_kernel("GPUTPCCFChargeMapFiller, fillIndexMap" "= TPCCLUSTERFINDER" LB single) o2_gpu_add_kernel("GPUTPCCFChargeMapFiller, fillFromDigits" "= TPCCLUSTERFINDER" LB single) diff --git a/GPU/Workflow/include/GPUWorkflow/GPUWorkflowSpec.h b/GPU/Workflow/include/GPUWorkflow/GPUWorkflowSpec.h index b218a21306a34..eda3b28c6cff6 100644 --- a/GPU/Workflow/include/GPUWorkflow/GPUWorkflowSpec.h +++ b/GPU/Workflow/include/GPUWorkflow/GPUWorkflowSpec.h @@ -233,6 +233,7 @@ class GPURecoWorkflowSpec : public o2::framework::Task bool mITSGeometryCreated = false; bool mTRDGeometryCreated = false; bool mPropagatorInstanceCreated = false; + int32_t mTPCCutAtTimeBin = -1; }; } // end namespace gpu diff --git a/GPU/Workflow/src/GPUWorkflowSpec.cxx b/GPU/Workflow/src/GPUWorkflowSpec.cxx index 94b1c3c2b8a7b..06942eab476c6 100644 --- a/GPU/Workflow/src/GPUWorkflowSpec.cxx +++ b/GPU/Workflow/src/GPUWorkflowSpec.cxx @@ -1013,9 +1013,8 @@ void GPURecoWorkflowSpec::doCalibUpdates(o2::framework::ProcessingContext& pc, c LOG(info) << "Updating solenoid field " << newCalibValues.solenoidField; } if (mAutoContinuousMaxTimeBin) { - mConfig->configGRP.grpContinuousMaxTimeBin = GPUO2InterfaceUtils::getTpcMaxTimeBinFromNHbf(mTFSettings->nHBFPerTF); newCalibValues.newContinuousMaxTimeBin = true; - newCalibValues.continuousMaxTimeBin = mConfig->configGRP.grpContinuousMaxTimeBin; + newCalibValues.continuousMaxTimeBin = mConfig->configGRP.grpContinuousMaxTimeBin = GPUO2InterfaceUtils::getTpcMaxTimeBinFromNHbf(mTFSettings->nHBFPerTF); LOG(info) << "Updating max time bin " << newCalibValues.continuousMaxTimeBin << " (" << mTFSettings->nHBFPerTF << " orbits)"; } @@ -1050,6 +1049,11 @@ void GPURecoWorkflowSpec::doCalibUpdates(o2::framework::ProcessingContext& pc, c if (mSpecConfig.runITSTracking) { needCalibUpdate = fetchCalibsCCDBITS(pc) || needCalibUpdate; } + if (mTPCCutAtTimeBin != mConfig->configGRP.tpcCutTimeBin) { + newCalibValues.newTPCTimeBinCut = true; + newCalibValues.tpcTimeBinCut = mConfig->configGRP.tpcCutTimeBin = mTPCCutAtTimeBin; + needCalibUpdate = true; + } if (needCalibUpdate) { LOG(info) << "Updating GPUReconstruction calibration objects"; mGPUReco->UpdateCalibration(newCalibObjects, newCalibValues); @@ -1095,8 +1099,13 @@ Inputs GPURecoWorkflowSpec::inputs() } else if (mSpecConfig.enableDoublePipeline == 1) { inputs.emplace_back("pipelineprepare", gDataOriginGPU, "PIPELINEPREPARE", 0, Lifetime::Timeframe); } + if (mSpecConfig.outputTracks || mSpecConfig.caClusterer) { + // calibration objects for TPC clusterization + inputs.emplace_back("tpcgain", gDataOriginTPC, "PADGAINFULL", 0, Lifetime::Condition, ccdbParamSpec(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalPadGainFull))); + inputs.emplace_back("tpcaltrosync", gDataOriginTPC, "ALTROSYNCSIGNAL", 0, Lifetime::Condition, ccdbParamSpec(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::AltroSyncSignal))); + } if (mSpecConfig.outputTracks) { - // loading calibration objects from the CCDB + // calibration objects for TPC tracking const auto mapSources = mSpecConfig.tpcDeadMapSources; if (mapSources != 0) { tpc::SourcesDeadMap sources((mapSources > -1) ? static_cast(mapSources) : tpc::SourcesDeadMap::All); @@ -1107,7 +1116,7 @@ Inputs GPURecoWorkflowSpec::inputs() inputs.emplace_back("tpcruninfo", gDataOriginTPC, "TPCRUNINFO", 0, Lifetime::Condition, ccdbParamSpec(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::ConfigRunInfo))); } } - inputs.emplace_back("tpcgain", gDataOriginTPC, "PADGAINFULL", 0, Lifetime::Condition, ccdbParamSpec(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalPadGainFull))); + inputs.emplace_back("tpcgainresidual", gDataOriginTPC, "PADGAINRESIDUAL", 0, Lifetime::Condition, ccdbParamSpec(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalPadGainResidual), {}, 1)); // time-dependent if (mSpecConfig.tpcUseMCTimeGain) { inputs.emplace_back("tpctimegain", gDataOriginTPC, "TIMEGAIN", 0, Lifetime::Condition, ccdbParamSpec(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalTimeGainMC), {}, 1)); // time-dependent @@ -1124,11 +1133,6 @@ Inputs GPURecoWorkflowSpec::inputs() if (mSpecConfig.decompressTPC) { inputs.emplace_back(InputSpec{"input", ConcreteDataTypeMatcher{gDataOriginTPC, mSpecConfig.decompressTPCFromROOT ? o2::header::DataDescription("COMPCLUSTERS") : o2::header::DataDescription("COMPCLUSTERSFLAT")}, Lifetime::Timeframe}); } else if (mSpecConfig.caClusterer) { - // if the output type are tracks, then the input spec for the gain map is already defined - if (!mSpecConfig.outputTracks) { - inputs.emplace_back("tpcgain", gDataOriginTPC, "PADGAINFULL", 0, Lifetime::Condition, ccdbParamSpec(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalPadGainFull))); - } - // We accept digits and MC labels also if we run on ZS Raw data, since they are needed for MC label propagation if ((!mSpecConfig.zsOnTheFly || mSpecConfig.processMC) && !mSpecConfig.zsDecoder) { inputs.emplace_back(InputSpec{"input", ConcreteDataTypeMatcher{gDataOriginTPC, "DIGITS"}, Lifetime::Timeframe}); diff --git a/GPU/Workflow/src/GPUWorkflowTPC.cxx b/GPU/Workflow/src/GPUWorkflowTPC.cxx index 97bf3aed26368..f895587b8b020 100644 --- a/GPU/Workflow/src/GPUWorkflowTPC.cxx +++ b/GPU/Workflow/src/GPUWorkflowTPC.cxx @@ -65,10 +65,7 @@ #include "SimulationDataFormat/MCCompLabel.h" #include "Algorithm/Parser.h" #include "DataFormatsGlobalTracking/RecoContainer.h" -#include "DataFormatsTRD/RecoInputContainer.h" -#include "TRDBase/Geometry.h" -#include "TRDBase/GeometryFlat.h" -#include "ITSBase/GeometryTGeo.h" +#include "DataFormatsTPC/AltroSyncSignal.h" #include "CommonUtils/VerbosityConfig.h" #include "CommonUtils/DebugStreamer.h" #include @@ -308,6 +305,10 @@ bool GPURecoWorkflowSpec::fetchCalibsCCDBTPC(ProcessingCon pc.inputs().get*>("tpcgain"); } + if (mSpecConfig.outputTracks || mSpecConfig.caClusterer) { + mTPCCutAtTimeBin = mConfParam->overrideTPCTimeBinCur > 0 ? mConfParam->overrideTPCTimeBinCur : pc.inputs().get("tpcaltrosync")->getTB2Cut(pc.services().get().tfCounter); + } + // these calibrations are only defined for the tracking if (mSpecConfig.outputTracks) { // update the calibration objects in case they changed in the CCDB diff --git a/Generators/CMakeLists.txt b/Generators/CMakeLists.txt index d60d185817c84..d909b3e604887 100644 --- a/Generators/CMakeLists.txt +++ b/Generators/CMakeLists.txt @@ -129,8 +129,15 @@ if(doBuildSimulation) COMPONENT_NAME Generator LABELS generator PUBLIC_LINK_LIBRARIES O2::Generators) + + o2_add_test(GeneratorPythia8Param NAME test_Generator_test_GeneratorPythia8Param + SOURCES test/test_GeneratorPythia8Param.cxx + COMPONENT_NAME Generator + LABELS generator + PUBLIC_LINK_LIBRARIES O2::Generators) endif() + o2_add_test_root_macro(share/external/tgenerator.C PUBLIC_LINK_LIBRARIES O2::Generators LABELS generators) diff --git a/Generators/include/Generators/GeneratorPythia8Param.h b/Generators/include/Generators/GeneratorPythia8Param.h index 165b1622239f5..612964fca73d9 100644 --- a/Generators/include/Generators/GeneratorPythia8Param.h +++ b/Generators/include/Generators/GeneratorPythia8Param.h @@ -24,28 +24,21 @@ namespace eventgen { /** - ** a parameter class/struct to keep the settings of - ** the Pythia8 event generator and - ** allow the user to modify them + ** a parameter class/struct to configure the settings of + ** the GeneratorPythia8 event generator **/ - -struct GeneratorPythia8Param : public o2::conf::ConfigurableParamHelper { +struct Pythia8GenConfig { std::string config = ""; std::string hooksFileName = ""; std::string hooksFuncName = ""; bool includePartonEvent = false; // whether to keep the event before hadronization std::string particleFilter = ""; // user particle filter int verbose = 0; // verbose control (if > 0 may show more info messages about what is going on) - O2ParamDef(GeneratorPythia8Param, "GeneratorPythia8"); }; -struct Pythia8GenConfig { - std::string config = ""; - std::string hooksFileName = ""; - std::string hooksFuncName = ""; - bool includePartonEvent = false; // whether to keep the event before hadronization - std::string particleFilter = ""; // user particle filter - int verbose = 0; // verbose control (if > 0 may show more info messages about what is going on) +// construct a configurable param singleton out of the Pythia8GenConfig struct +struct GeneratorPythia8Param : public o2::conf::ConfigurableParamPromoter { + O2ParamDef(GeneratorPythia8Param, "GeneratorPythia8"); }; } // end namespace eventgen diff --git a/Generators/src/GeneratorPythia8.cxx b/Generators/src/GeneratorPythia8.cxx index 8c9b4fcffdff2..fef2c4d2e9a1c 100644 --- a/Generators/src/GeneratorPythia8.cxx +++ b/Generators/src/GeneratorPythia8.cxx @@ -45,26 +45,11 @@ namespace eventgen /*****************************************************************/ /*****************************************************************/ -GeneratorPythia8::GeneratorPythia8() : Generator("ALICEo2", "ALICEo2 Pythia8 Generator") +// the default construct uses the GeneratorPythia8Param singleton to extract a config and delegates +// to the proper constructor +GeneratorPythia8::GeneratorPythia8() : GeneratorPythia8(GeneratorPythia8Param::Instance().detach()) { - /** default constructor **/ - - mInterface = reinterpret_cast(&mPythia); - mInterfaceName = "pythia8"; - - auto& param = GeneratorPythia8Param::Instance(); - LOG(info) << "Default Instance \'Pythia8\' generator with following parameters"; - LOG(info) << param; - - // convert the outside singleton config to the internally used one - o2::eventgen::Pythia8GenConfig config{param.config, - param.hooksFileName, param.hooksFuncName, param.includePartonEvent, param.particleFilter, param.verbose}; - mGenConfig = config; - - setConfig(config.config); - setHooksFileName(config.hooksFileName); - setHooksFuncName(config.hooksFuncName); - // TODO: use constructor delegation to other interface + LOG(info) << "GeneratorPythia8 constructed from GeneratorPythia8Param ConfigurableParam"; } /*****************************************************************/ diff --git a/Generators/src/GeneratorPythia8Param.cxx b/Generators/src/GeneratorPythia8Param.cxx index 984680e46ad01..6b477beb16ba9 100644 --- a/Generators/src/GeneratorPythia8Param.cxx +++ b/Generators/src/GeneratorPythia8Param.cxx @@ -12,4 +12,4 @@ /// \author R+Preghenella - January 2020 #include "Generators/GeneratorPythia8Param.h" -O2ParamImpl(o2::eventgen::GeneratorPythia8Param); +O2ParamImpl(o2::eventgen::GeneratorPythia8Param); \ No newline at end of file diff --git a/Generators/src/GeneratorsLinkDef.h b/Generators/src/GeneratorsLinkDef.h index 41e14b02f18b9..fe219c6f5476c 100644 --- a/Generators/src/GeneratorsLinkDef.h +++ b/Generators/src/GeneratorsLinkDef.h @@ -48,8 +48,9 @@ #pragma link C++ class o2::eventgen::GeneratorPythia8Param + ; #pragma link C++ class o2::eventgen::Pythia8GenConfig + ; #pragma link C++ class o2::eventgen::DecayerPythia8Param + ; -#pragma link C++ class o2::conf::ConfigurableParamHelper < o2::eventgen::GeneratorPythia8Param> + ; +#pragma link C++ class o2::conf::ConfigurableParamPromoter < o2::eventgen::GeneratorPythia8Param, o2::eventgen::Pythia8GenConfig> + ; #pragma link C++ class o2::conf::ConfigurableParamHelper < o2::eventgen::DecayerPythia8Param> + ; + #pragma link C++ class o2::eventgen::GeneratorFactory + ; #endif #if defined(GENERATORS_WITH_PYTHIA8) && defined(GENERATORS_WITH_HEPMC3) diff --git a/Generators/test/test_GeneratorPythia8Param.cxx b/Generators/test/test_GeneratorPythia8Param.cxx new file mode 100644 index 0000000000000..c735487ea293c --- /dev/null +++ b/Generators/test/test_GeneratorPythia8Param.cxx @@ -0,0 +1,80 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#define BOOST_TEST_MODULE Test GeneratorPythia8Param class +#define BOOST_TEST_MAIN +#define BOOST_TEST_DYN_LINK +#include +#include +#include +#include +#include "CCDB/BasicCCDBManager.h" + +// Tests various aspects of the +// ConfigurableParamPromoter class, which is used to promote +// Pythia8GenConfig to a configurable param +BOOST_AUTO_TEST_CASE(pythia8_Pythia8GenConfig) +{ + o2::conf::ConfigurableParam::updateFromString( + "GeneratorPythia8.config=Foo;GeneratorPythia8.includePartonEvent=true"); + + using o2::eventgen::GeneratorPythia8Param; + + BOOST_CHECK_EQUAL(GeneratorPythia8Param::Instance().config, std::string("Foo")); + BOOST_CHECK_EQUAL(GeneratorPythia8Param::Instance().includePartonEvent, true); + + BOOST_CHECK_EQUAL(GeneratorPythia8Param::Instance().includePartonEvent, o2::conf::ConfigurableParam::getValueAs("GeneratorPythia8.includePartonEvent")); + // setValue - getValue + o2::conf::ConfigurableParam::setValue("GeneratorPythia8.config", "Baz"); + BOOST_CHECK_EQUAL(o2::conf::ConfigurableParam::getValueAs("GeneratorPythia8.config"), std::string("Baz")); + BOOST_CHECK_EQUAL(GeneratorPythia8Param::Instance().config, std::string("Baz")); + + // member provenance + BOOST_CHECK_EQUAL(GeneratorPythia8Param::Instance().getMemberProvenance("config"), o2::conf::ConfigurableParam::EParamProvenance::kRT); + BOOST_CHECK_EQUAL(GeneratorPythia8Param::Instance().getMemberProvenance("verbose"), o2::conf::ConfigurableParam::EParamProvenance::kCODE); + + // config detach + auto config_copy = GeneratorPythia8Param::Instance().detach(); + BOOST_CHECK_EQUAL(config_copy.config, std::string("Baz")); + BOOST_CHECK_EQUAL(config_copy.includePartonEvent, true); + + // file IO + TFile tmp_file("GeneratorParamConfig_tmp.root", "RECREATE"); + + GeneratorPythia8Param::Instance().serializeTo(&tmp_file); + // modify the instance to some intermediate fluent value + o2::conf::ConfigurableParam::setValue("GeneratorPythia8.includePartonEvent", "0"); + BOOST_CHECK_EQUAL(config_copy.includePartonEvent, true); + BOOST_CHECK_EQUAL(GeneratorPythia8Param::Instance().includePartonEvent, false); + tmp_file.Close(); + + // read back + TFile tmp_file2("GeneratorParamConfig_tmp.root", "READ"); + const_cast(GeneratorPythia8Param::Instance()).initFrom(&tmp_file2); + BOOST_CHECK_EQUAL(GeneratorPythia8Param::Instance().includePartonEvent, true); + tmp_file2.Close(); + + // CCDB IO + std::string ccdbUrl = "http://ccdb-test.cern.ch:8080"; + bool hostReachable = false; + o2::ccdb::CcdbApi api; + api.init(ccdbUrl); + std::string pathA = "/Generators/UnitTest/Pythia8/GeneratorPythia8Param"; + std::map md; + long start = 1000, stop = 2000; + api.storeAsTFileAny(&GeneratorPythia8Param::Instance(), pathA, md, start, stop); + + // modify the instance to some intermediate fluent value + o2::conf::ConfigurableParam::setValue("GeneratorPythia8.includePartonEvent", "0"); + + auto returnedobj = api.retrieveFromTFileAny(pathA, md, (start + stop) / 2); + GeneratorPythia8Param::Instance().printKeyValues(); +}; diff --git a/Steer/src/O2MCApplication.cxx b/Steer/src/O2MCApplication.cxx index 96cc2f2e969db..02d332b0c0641 100644 --- a/Steer/src/O2MCApplication.cxx +++ b/Steer/src/O2MCApplication.cxx @@ -264,29 +264,19 @@ void addSpecialParticles() LOG(info) << "Adding custom particles to VMC"; //Hypertriton - TVirtualMC::GetMC()->DefineParticle(1010010030, "HyperTriton", kPTHadron, 2.99131, 1.0, 2.632e-10, "Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 3, kFALSE); + TVirtualMC::GetMC()->DefineParticle(1010010030, "HyperTriton", kPTHadron, 2.991134, 1.0, 2.632e-10, "Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 3, kFALSE); //Anti-Hypertriton - TVirtualMC::GetMC()->DefineParticle(-1010010030, "AntiHyperTriton", kPTHadron, 2.99131, 1.0, 2.632e-10, "Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 3, kFALSE); + TVirtualMC::GetMC()->DefineParticle(-1010010030, "AntiHyperTriton", kPTHadron, 2.991134, 1.0, 2.632e-10, "Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 3, kFALSE); //Hyper hydrogen 4 ground state - TVirtualMC::GetMC()->DefineParticle(1010010040, "Hyperhydrog4", kPTHadron, 3.9226, 1.0, 2.632e-10, "Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 4, kFALSE); + TVirtualMC::GetMC()->DefineParticle(1010010040, "Hyperhydrog4", kPTHadron, 3.922434, 1.0, 2.08e-10, "Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 4, kFALSE); //Anti-Hyper hydrogen 4 ground state - TVirtualMC::GetMC()->DefineParticle(-1010010040, "AntiHyperhydrog4", kPTHadron, 3.9226, 1.0, 2.632e-10, "Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 4, kFALSE); - - //Hyper hydrogen 4 excited state - TVirtualMC::GetMC()->DefineParticle(1010010041, "Hyperhydrog4*", kPTHadron, 3.9237, 1.0, 2.632e-10, "Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 4, kFALSE); - //Anti-Hyper hydrogen 4 excited state - TVirtualMC::GetMC()->DefineParticle(-1010010041, "AntiHyperhydrog4*", kPTHadron, 3.9237, 1.0, 2.632e-10, "Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 4, kFALSE); + TVirtualMC::GetMC()->DefineParticle(-1010010040, "AntiHyperhydrog4", kPTHadron, 3.922434, 1.0, 2.08e-10, "Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 4, kFALSE); //Hyper helium 4 ground state - TVirtualMC::GetMC()->DefineParticle(1010020040, "Hyperhelium4", kPTHadron, 3.9217, 2.0, 2.632e-10, "Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 4, kFALSE); + TVirtualMC::GetMC()->DefineParticle(1010020040, "Hyperhelium4", kPTHadron, 3.921728, 2.0, 2.50e-10, "Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 4, kFALSE); //Anti-Hyper helium 4 ground state - TVirtualMC::GetMC()->DefineParticle(-1010020040, "AntiHyperhelium4", kPTHadron, 3.9217, 2.0, 2.632e-10, "Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 4, kFALSE); - - //Hyper helium 4 excited state - TVirtualMC::GetMC()->DefineParticle(1010020041, "Hyperhelium4*", kPTHadron, 3.9231, 2.0, 2.632e-10, "Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 4, kFALSE); - //Anti-Hyper helium 4 excited state - TVirtualMC::GetMC()->DefineParticle(-1010020041, "AntiHyperhelium4*", kPTHadron, 3.9231, 2.0, 2.632e-10, "Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 4, kFALSE); + TVirtualMC::GetMC()->DefineParticle(-1010020040, "AntiHyperhelium4", kPTHadron, 3.921728, 2.0, 2.50e-10, "Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 4, kFALSE); // Lithium 4 ground state TVirtualMC::GetMC()->DefineParticle(1000030040, "Lithium4", kPTHadron, 3.7513, 3.0, 9.1e-23, "Ion", 0.003, 0, 1, 0, 0, 0, 0, 0, 4, kFALSE); @@ -294,15 +284,24 @@ void addSpecialParticles() TVirtualMC::GetMC()->DefineParticle(-1000030040, "AntiLithium4", kPTHadron, 3.7513, 3.0, 9.1e-23, "Ion", 0.003, 0, 1, 0, 0, 0, 0, 0, 4, kFALSE); //Hyper helium 5 - TVirtualMC::GetMC()->DefineParticle(1010020050, "Hyperhelium5", kPTHadron, 4.841, 2.0, 2.632e-10, "Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 5, kFALSE); + TVirtualMC::GetMC()->DefineParticle(1010020050, "Hyperhelium5", kPTHadron, 4.839961, 2.0, 2.74e-10, "Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 5, kFALSE); //Anti-Hyper helium 5 - TVirtualMC::GetMC()->DefineParticle(-1010020050, "AntiHyperhelium5", kPTHadron, 4.841, 2.0, 2.632e-10, "Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 5, kFALSE); + TVirtualMC::GetMC()->DefineParticle(-1010020050, "AntiHyperhelium5", kPTHadron, 4.839961, 2.0, 2.74e-10, "Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 5, kFALSE); //Double Hyper hydrogen 4 TVirtualMC::GetMC()->DefineParticle(1020010040, "DoubleHyperhydrogen4", kPTHadron, 4.106, 1.0, 2.632e-10, "Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 4, kFALSE); //Double Anti-Hyper hydrogen 4 TVirtualMC::GetMC()->DefineParticle(-1020010040, "DoubleAntiHyperhydrogen4", kPTHadron, 4.106, 1.0, 2.632e-10, "Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 4, kFALSE); + // 4Xi(-)H + TVirtualMC::GetMC()->DefineParticle(1120010040, "4XiH", kPTHadron, 4.128, 1.0, 1.639e-10, "Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 4, kFALSE); + // Anti-4Xi(-)H + TVirtualMC::GetMC()->DefineParticle(-1120010040, "Anti4XiH", kPTHadron, 4.128, 1.0, 1.639e-10, "Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 4, kFALSE); + // 4Xi(-)He + TVirtualMC::GetMC()->DefineParticle(1120020040, "4XiHe", kPTHadron, 4.128, 1.0, 1.639e-10, "Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 4, kFALSE); + // Anti-4Xi(-)He + TVirtualMC::GetMC()->DefineParticle(-1120020040, "Anti4XiHe", kPTHadron, 4.128, 1.0, 1.639e-10, "Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 4, kFALSE); + // Hyper helium 4 sigma TVirtualMC::GetMC()->DefineParticle(1110020040, "Hyperhelium4sigma", kPTHadron, 3.995, 2.0, 8.018e-11, "Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 4, kFALSE); // Anti-Hyper helium 4 sigma @@ -586,8 +585,6 @@ void addSpecialParticles() mode3[1][2] = -211; // negative pion TVirtualMC::GetMC()->SetDecayMode(1010010040, bratio3, mode3); - //Decay for the excited state (after em transition) - TVirtualMC::GetMC()->SetDecayMode(1010010041, bratio3, mode3); // Define the 2- and 3-body phase space decay for the Hyper Hydrogen 4 Int_t amode3[6][3]; @@ -608,8 +605,6 @@ void addSpecialParticles() amode3[1][2] = 211; // positive pion TVirtualMC::GetMC()->SetDecayMode(-1010010040, abratio3, amode3); - //Decay for the excited state (after em transition) - TVirtualMC::GetMC()->SetDecayMode(-1010010041, abratio3, amode3); // Define the 3-body phase space decay for the Hyper Helium 4 Int_t mode4[6][3]; @@ -621,14 +616,16 @@ void addSpecialParticles() mode4[kz][1] = 0; mode4[kz][2] = 0; } - bratio4[0] = 100.; + bratio4[0] = 50.; mode4[0][0] = 1000020030; // Helium3 mode4[0][1] = -211; // negative pion mode4[0][2] = 2212; // proton + bratio4[1] = 50.; + mode4[1][0] = 1000030040; // lithium-4 + mode4[1][1] = -211; // negative pion + TVirtualMC::GetMC()->SetDecayMode(1010020040, bratio4, mode4); - //Decay for the excited state (after em transition) - TVirtualMC::GetMC()->SetDecayMode(1010020041, bratio4, mode4); // Define the 2-body phase space decay for the Anti-Hyper Helium 4 Int_t amode4[6][3]; @@ -640,14 +637,16 @@ void addSpecialParticles() amode4[kz][1] = 0; amode4[kz][2] = 0; } - abratio4[0] = 100.; + abratio4[0] = 50.; amode4[0][0] = -1000020030; // anti-Helium 3 amode4[0][1] = 211; // positive pion amode4[0][2] = -2212; // anti proton + abratio4[1] = 50.; + amode4[1][0] = -1000030040; // antilithium-4 + amode4[1][1] = 211; // positive pion + TVirtualMC::GetMC()->SetDecayMode(-1010020040, abratio4, amode4); - //Decay for the excited state (after em transition) - TVirtualMC::GetMC()->SetDecayMode(-1010020041, abratio4, amode4); // Define the 2-body phase space decay for the Lithium 4 Int_t model4[6][3]; @@ -733,10 +732,15 @@ void addSpecialParticles() mode42[kz][1] = 0; mode42[kz][2] = 0; } - bratio42[0] = 100.; + bratio42[0] = 50.; mode42[0][0] = 1010020040; // Hyper-Helium4 mode42[0][1] = -211; // negative pion + bratio42[1] = 50.; + mode42[1][0] = 1010010030; // Hypertriton + mode42[1][1] = 2212; // proton + mode42[1][2] = -211; // negative pion + TVirtualMC::GetMC()->SetDecayMode(1020010040, bratio42, mode42); // Define the 2-body phase space decay for the Anti Double Hyper Hydrogen 4 @@ -749,12 +753,121 @@ void addSpecialParticles() amode42[kz][1] = 0; amode42[kz][2] = 0; } - abratio42[0] = 100.; + abratio42[0] = 50.; amode42[0][0] = -1010020040; // anti-Hyper-Helium 4 amode42[0][1] = 211; // positive pion + abratio42[1] = 50.; + amode42[1][0] = -1010010030; // anti-Hypertriton + amode42[1][1] = -2212; // antiproton + amode42[1][2] = 211; // positive pion + TVirtualMC::GetMC()->SetDecayMode(-1020010040, abratio42, amode42); + // Define the decay for the 4Xi(-)He + Int_t mode4XiHe[6][3]; + Float_t bratio4XiHe[6]; + + for (Int_t kz = 0; kz < 6; kz++) { + bratio4XiHe[kz] = 0.; + mode4XiHe[kz][0] = 0; + mode4XiHe[kz][1] = 0; + mode4XiHe[kz][2] = 0; + } + bratio4XiHe[0] = 33.; + mode4XiHe[0][0] = 1010020040; // HyperHelium4 + mode4XiHe[0][1] = -211; // negative pion + + bratio4XiHe[1] = 33.; + mode4XiHe[1][0] = 3122; // lambda + mode4XiHe[1][1] = 1000020030; // helium-3 + mode4XiHe[1][2] = -211; // negative pion + + bratio4XiHe[2] = 33.; + mode4XiHe[2][0] = 1000030040; // lithium-4 + mode4XiHe[2][1] = -211; // negative pion + mode4XiHe[2][2] = -211; // negative pion + + TVirtualMC::GetMC()->SetDecayMode(1120020040, bratio4XiHe, mode4XiHe); + + // Define the decay for the Anti-4Xi(-)He + Int_t amode4XiHe[6][3]; + Float_t abratio4XiHe[6]; + + for (Int_t kz = 0; kz < 6; kz++) { + abratio4XiHe[kz] = 0.; + amode4XiHe[kz][0] = 0; + amode4XiHe[kz][1] = 0; + amode4XiHe[kz][2] = 0; + } + abratio4XiHe[0] = 33.; + amode4XiHe[0][0] = -1010020040; // antiHyperHelium-4 + amode4XiHe[0][1] = 211; // positive pion + + abratio4XiHe[1] = 33.; + amode4XiHe[1][0] = -3122; // antilambda + amode4XiHe[1][1] = -1000020030; // antihelium-3 + amode4XiHe[1][2] = 211; // positive pion + + abratio4XiHe[2] = 33.; + amode4XiHe[2][0] = -1000030040; // antilithium-4 + amode4XiHe[2][1] = 211; // positive pion + amode4XiHe[2][2] = 211; // positive pion + + TVirtualMC::GetMC()->SetDecayMode(-1120020040, abratio4XiHe, amode4XiHe); + + // Define the decay for the 4Xi(-)H + Int_t mode4XiH[6][3]; + Float_t bratio4XiH[6]; + + for (Int_t kz = 0; kz < 6; kz++) { + bratio4XiH[kz] = 0.; + mode4XiH[kz][0] = 0; + mode4XiH[kz][1] = 0; + mode4XiH[kz][2] = 0; + } + bratio4XiH[0] = 33.; + mode4XiH[0][0] = 1010010040; // HyperHydrogen4 + mode4XiH[0][1] = -211; // negative pion + + bratio4XiH[1] = 33.; + mode4XiH[1][0] = 3122; // lambda + mode4XiH[1][1] = 1000010030; // triton + mode4XiH[1][2] = -211; // negative pion + + bratio4XiH[2] = 33.; + mode4XiH[2][0] = 1000020040; // alpha + mode4XiH[2][1] = -211; // negative pion + mode4XiH[2][2] = -211; // negative pion + + TVirtualMC::GetMC()->SetDecayMode(1120010040, bratio4XiH, mode4XiH); + + // Define the decay for the Anti-4Xi(-)H + Int_t amode4XiH[6][3]; + Float_t abratio4XiH[6]; + + for (Int_t kz = 0; kz < 6; kz++) { + abratio4XiH[kz] = 0.; + amode4XiH[kz][0] = 0; + amode4XiH[kz][1] = 0; + amode4XiH[kz][2] = 0; + } + abratio4XiH[0] = 33.; + amode4XiH[0][0] = -1010010040; // antiHyperHydrogen-4 + amode4XiH[0][1] = 211; // positive pion + + abratio4XiH[1] = 33.; + amode4XiH[1][0] = -3122; // antilambda + amode4XiH[1][1] = -1000010030; // antitriton + amode4XiH[1][2] = 211; // positive pion + + abratio4XiH[2] = 33.; + amode4XiH[2][0] = -1000020040; // antialpha + amode4XiH[2][1] = 211; // positive pion + amode4XiH[2][2] = 211; // positive pion + + TVirtualMC::GetMC()->SetDecayMode(-1120010040, abratio4XiH, amode4XiH); + // Define the 2- and 3-body phase space decay for the Hyper Helium 4 sigma Int_t mode4s[6][3]; Float_t bratio4s[6]; diff --git a/run/CMakeLists.txt b/run/CMakeLists.txt index f21ecafb0528a..662716901ed0a 100644 --- a/run/CMakeLists.txt +++ b/run/CMakeLists.txt @@ -193,6 +193,8 @@ o2_add_test_command(NAME o2sim_G4 2 --skipModules MFT ZDC + --seed + 15946057944514955802 --configKeyValues "align-geom.mDetectors=none" ENVIRONMENT "${SIMENV}" @@ -255,6 +257,8 @@ o2_add_test_command(NAME o2sim_G3 pythia8pp --chunkSize 10 + --seed + 15946057944514955802 --configKeyValues "align-geom.mDetectors=none" LABELS g3 sim long diff --git a/run/checkStack.cxx b/run/checkStack.cxx index 4470ea463fd98..36a9da7a62c13 100644 --- a/run/checkStack.cxx +++ b/run/checkStack.cxx @@ -142,10 +142,10 @@ int main(int argc, char** argv) bool havereferences = trackrefs->size(); if (havereferences) { for (auto& trackID : trackidsinTPC) { - auto trackrefs = mcreader.getTrackRefs(eventID, trackID); - assert(trackrefs.size() > 0); - LOG(debug) << " Track " << trackID << " has " << trackrefs.size() << " TrackRefs"; - for (auto& ref : trackrefs) { + auto tpc_trackrefs = mcreader.getTrackRefs(eventID, trackID); + LOG(debug) << " Track " << trackID << " has " << tpc_trackrefs.size() << " TrackRefs"; + assert(tpc_trackrefs.size() > 0); + for (auto& ref : tpc_trackrefs) { assert(ref.getTrackID() == trackID); } }