diff --git a/Algorithm/include/Algorithm/PageParser.h b/Algorithm/include/Algorithm/PageParser.h index e382fc318352e..3ca01d87bcba3 100644 --- a/Algorithm/include/Algorithm/PageParser.h +++ b/Algorithm/include/Algorithm/PageParser.h @@ -255,12 +255,12 @@ class PageParser return mElement; } // comparison - bool operator==(const SelfType& rh) + bool operator==(const SelfType& rh) const { return mPosition == rh.mPosition; } // comparison - bool operator!=(const SelfType& rh) + bool operator!=(const SelfType& rh) const { return mPosition != rh.mPosition; } diff --git a/CCDB/include/CCDB/CCDBDownloader.h b/CCDB/include/CCDB/CCDBDownloader.h index f2706bd26e580..4389dc54c4ce3 100644 --- a/CCDB/include/CCDB/CCDBDownloader.h +++ b/CCDB/include/CCDB/CCDBDownloader.h @@ -51,6 +51,7 @@ typedef struct DownloaderRequestData { long timestamp; HeaderObjectPair_t hoPair; std::map* headers; + std::string userAgent; std::function localContentCallback; bool errorflag = false; @@ -208,6 +209,11 @@ class CCDBDownloader */ void runLoop(bool noWait); + /** + * Returns a message describing the transfer an it's result. + */ + std::string prepareLogMessage(std::string host_url, std::string userAgent, const std::string& path, long ts, const std::map* headers, long httpCode) const; + /** * Leaves only the protocol and host part of the url, discrading path and metadata. */ diff --git a/CCDB/src/CCDBDownloader.cxx b/CCDB/src/CCDBDownloader.cxx index b08fe176505ff..c3cde9cd8f469 100644 --- a/CCDB/src/CCDBDownloader.cxx +++ b/CCDB/src/CCDBDownloader.cxx @@ -493,21 +493,25 @@ void CCDBDownloader::transferFinished(CURL* easy_handle, CURLcode curlCode) curl_easy_getinfo(easy_handle, CURLINFO_EFFECTIVE_URL, &url); LOG(debug) << "Transfer for " << url << " finished with code " << httpCode << "\n"; + std::string currentHost = requestData->hosts[performData->hostInd]; + std::string loggingMessage = prepareLogMessage(currentHost, requestData->userAgent, requestData->path, requestData->timestamp, requestData->headers, httpCode); + // Get alternative locations for the same host auto locations = getLocations(&(requestData->hoPair.header)); // React to received http code - if (404 == httpCode) { - LOG(error) << "Requested resource does not exist: " << url; - } else if (304 == httpCode) { - LOGP(debug, "Object exists but I am not serving it since it's already in your possession"); - contentRetrieved = true; - } else if (300 <= httpCode && httpCode < 400 && performData->locInd < locations.size()) { - followRedirect(performData, easy_handle, locations, rescheduled, contentRetrieved); - } else if (200 <= httpCode && httpCode < 300) { - contentRetrieved = true; + if (200 <= httpCode && httpCode < 400) { + LOG(debug) << loggingMessage; + if (304 == httpCode) { + LOGP(debug, "Object exists but I am not serving it since it's already in your possession"); + contentRetrieved = true; + } else if (300 <= httpCode && httpCode < 400 && performData->locInd < locations.size()) { + followRedirect(performData, easy_handle, locations, rescheduled, contentRetrieved); + } else if (200 <= httpCode && httpCode < 300) { + contentRetrieved = true; + } } else { - LOG(error) << "Error in fetching object " << url << ", curl response code:" << httpCode; + LOG(error) << loggingMessage; } // Check if content was retrieved, or scheduled to be retrieved @@ -694,4 +698,21 @@ void CCDBDownloader::asynchSchedule(CURL* handle, size_t* requestCounter) // return codeVector; } +std::string CCDBDownloader::prepareLogMessage(std::string host_url, std::string userAgent, const std::string& path, long ts, const std::map* headers, long httpCode) const +{ + std::string upath{path}; + if (headers) { + auto ent = headers->find("Valid-From"); + if (ent != headers->end()) { + upath += "/" + ent->second; + } + ent = headers->find("ETag"); + if (ent != headers->end()) { + upath += "/" + ent->second; + } + } + upath.erase(remove(upath.begin(), upath.end(), '\"'), upath.end()); + return fmt::format("CcdbDownloader finished transfer {}{}{} for {} (agent_id: {}) with http code: {}", host_url, (host_url.back() == '/') ? "" : "/", upath, (ts < 0) ? getCurrentTimestamp() : ts, userAgent, httpCode); +} + } // namespace o2 diff --git a/CCDB/src/CcdbApi.cxx b/CCDB/src/CcdbApi.cxx index 6bcc7be9f2988..cc10effec1c6e 100644 --- a/CCDB/src/CcdbApi.cxx +++ b/CCDB/src/CcdbApi.cxx @@ -1533,6 +1533,7 @@ void CcdbApi::scheduleDownload(RequestContext& requestContext, size_t* requestCo data->path = requestContext.path; data->timestamp = requestContext.timestamp; data->localContentCallback = localContentCallback; + data->userAgent = mUniqueAgentID; curl_easy_setopt(curl_handle, CURLOPT_URL, fullUrl.c_str()); initCurlOptionsForRetrieve(curl_handle, (void*)(&data->hoPair), writeCallback, false); @@ -1671,8 +1672,6 @@ void CcdbApi::vectoredLoadFileToMemory(std::vector& requestConte // navigateSourcesAndLoadFile either retrieves file from snapshot immediately, or schedules it to be downloaded when mDownloader->runLoop is ran at a later time auto& requestContext = requestContexts.at(i); navigateSourcesAndLoadFile(requestContext, fromSnapshots.at(i), &requestCounter); - logReading(requestContext.path, requestContext.timestamp, &requestContext.headers, - fmt::format("{}{}", requestContext.considerSnapshot ? "load to memory" : "retrieve", fromSnapshots.at(i) ? " from snapshot" : "")); } // Download the rest @@ -1683,12 +1682,12 @@ void CcdbApi::vectoredLoadFileToMemory(std::vector& requestConte // Save snapshots for (int i = 0; i < requestContexts.size(); i++) { auto& requestContext = requestContexts.at(i); + logReading(requestContext.path, requestContext.timestamp, &requestContext.headers, + fmt::format("{}{}", requestContext.considerSnapshot ? "load to memory" : "retrieve", fromSnapshots.at(i) ? " from snapshot" : "")); if (!requestContext.dest.empty()) { if (requestContext.considerSnapshot && fromSnapshots.at(i) != 2) { saveSnapshot(requestContext); } - } else { - LOG(warning) << "Did not receive content for " << requestContext.path << "\n"; // Temporarily demoted to warning, since it floods the infologger } } } diff --git a/Common/Constants/include/CommonConstants/PhysicsConstants.h b/Common/Constants/include/CommonConstants/PhysicsConstants.h index 58a7300dc4b98..56357ab5ea5d6 100644 --- a/Common/Constants/include/CommonConstants/PhysicsConstants.h +++ b/Common/Constants/include/CommonConstants/PhysicsConstants.h @@ -12,40 +12,167 @@ /// \file PhysicsConstants.h /// \brief Header to collect physics constants /// \author ruben.shahoyan@cern.ch +/// \author Vít Kučera , Inha University +/// \note Use the make_pdg_header.py script to generate the enums and mass declarations. #ifndef ALICEO2_PHYSICSCONSTANTS_H_ #define ALICEO2_PHYSICSCONSTANTS_H_ -namespace o2 -{ -namespace constants -{ -namespace physics +namespace o2::constants::physics { // particles masses -constexpr float MassPhoton = 0.0; -constexpr float MassElectron = 0.000511; -constexpr float MassMuon = 0.105658; -constexpr float MassPionCharged = 0.139570; -constexpr float MassPionNeutral = 0.134976; -constexpr float MassKaonCharged = 0.493677; -constexpr float MassKaonNeutral = 0.497648; -constexpr float MassProton = 0.938272; -constexpr float MassLambda = 1.115683; -constexpr float MassDeuteron = 1.8756129; -constexpr float MassTriton = 2.8089211; -constexpr float MassHelium3 = 2.8083916; -constexpr float MassAlpha = 3.7273794; -constexpr float MassHyperTriton = 2.992; -constexpr float MassHyperhydrog4 = 3.931; -constexpr float MassHyperhelium4 = 3.9218; -constexpr float MassXiMinus = 1.32171; -constexpr float MassOmegaMinus = 1.67245; + +// BEGINNING OF THE GENERATED BLOCK. +// DO NOT EDIT THIS BLOCK DIRECTLY! +// It has been generated by the make_pdg_header.py script. +// For modifications, edit the script and generate this block again. + +/// \brief Declarations of named PDG codes of particles missing in ROOT PDG_t +/// \note Follow kCamelCase naming convention +/// \link https://root.cern/doc/master/TPDGCode_8h.html +enum Pdg { + kB0 = 511, + kB0Bar = -511, + kBPlus = 521, + kBS = 531, + kBSBar = -531, + kD0 = 421, + kD0Bar = -421, + kDMinus = -411, + kDPlus = 411, + kDS = 431, + kDSBar = -431, + kDStar = 413, + kChiC1 = 20443, + kJPsi = 443, + kLambdaB0 = 5122, + kLambdaCPlus = 4122, + kOmegaC0 = 4332, + kPhi = 333, + kSigmaC0 = 4112, + kSigmaCPlusPlus = 4222, + kX3872 = 9920443, + kXi0 = 3322, + kXiB0 = 5232, + kXiCCPlusPlus = 4422, + kXiCPlus = 4232, + kXiCZero = 4132, + kDeuteron = 1000010020, + kTriton = 1000010030, + kHelium3 = 1000020030, + kAlpha = 1000020040, + kHyperTriton = 1010010030, + kHyperHydrogen4 = 1010010040, + kHyperHelium4 = 1010020040 +}; + +/// \brief Declarations of masses for additional particles +constexpr double MassB0 = 5.27953; +constexpr double MassB0Bar = 5.27953; +constexpr double MassBPlus = 5.27915; +constexpr double MassBS = 5.3663; +constexpr double MassBSBar = 5.3663; +constexpr double MassD0 = 1.86484; +constexpr double MassD0Bar = 1.86484; +constexpr double MassDMinus = 1.86962; +constexpr double MassDPlus = 1.86962; +constexpr double MassDS = 1.9685; +constexpr double MassDSBar = 1.9685; +constexpr double MassDStar = 2.01027; +constexpr double MassChiC1 = 3.51066; +constexpr double MassJPsi = 3.096916; +constexpr double MassLambdaB0 = 5.6202; +constexpr double MassLambdaCPlus = 2.28646; +constexpr double MassOmegaC0 = 2.6975; +constexpr double MassPhi = 1.019455; +constexpr double MassSigmaC0 = 2.45376; +constexpr double MassSigmaCPlusPlus = 2.45402; +constexpr double MassX3872 = 3.87165; +constexpr double MassXi0 = 1.31486; +constexpr double MassXiB0 = 5.7924; +constexpr double MassXiCCPlusPlus = 3.59798; +constexpr double MassXiCPlus = 2.4679; +constexpr double MassXiCZero = 2.471; +constexpr double MassDeuteron = 1.87561294257; +constexpr double MassTriton = 2.80892113298; +constexpr double MassHelium3 = 2.80839160743; +constexpr double MassAlpha = 3.7273794066; +constexpr double MassHyperTriton = 2.99131; +constexpr double MassHyperHydrogen4 = 3.9226; +constexpr double MassHyperHelium4 = 3.9217; + +/// \brief Declarations of masses for particles in ROOT PDG_t +constexpr double MassDown = 0.0048; +constexpr double MassDownBar = 0.0048; +constexpr double MassUp = 0.0024; +constexpr double MassUpBar = 0.0024; +constexpr double MassStrange = 0.104; +constexpr double MassStrangeBar = 0.104; +constexpr double MassCharm = 1.27; +constexpr double MassCharmBar = 1.27; +constexpr double MassBottom = 4.68; +constexpr double MassBottomBar = 4.68; +constexpr double MassTop = 171.2; +constexpr double MassTopBar = 171.2; +constexpr double MassGluon = 0.0; +constexpr double MassElectron = 0.00051099891; +constexpr double MassPositron = 0.00051099891; +constexpr double MassNuE = 0.0; +constexpr double MassNuEBar = 0.0; +constexpr double MassMuonMinus = 0.105658; +constexpr double MassMuonPlus = 0.105658; +constexpr double MassNuMu = 0.0; +constexpr double MassNuMuBar = 0.0; +constexpr double MassTauMinus = 1.77684; +constexpr double MassTauPlus = 1.77684; +constexpr double MassNuTau = 0.0; +constexpr double MassNuTauBar = 0.0; +constexpr double MassGamma = 0.0; +constexpr double MassZ0 = 91.187; +constexpr double MassWPlus = 80.398; +constexpr double MassWMinus = 80.398; +constexpr double MassPi0 = 0.134977; +constexpr double MassK0Long = 0.497614; +constexpr double MassPiPlus = 0.13957; +constexpr double MassPiMinus = 0.13957; +constexpr double MassProton = 0.938272; +constexpr double MassProtonBar = 0.938272; +constexpr double MassNeutron = 0.939565; +constexpr double MassNeutronBar = 0.939565; +constexpr double MassK0Short = 0.497614; +constexpr double MassK0 = 0.497614; +constexpr double MassK0Bar = 0.497614; +constexpr double MassKPlus = 0.493677; +constexpr double MassKMinus = 0.493677; +constexpr double MassLambda0 = 1.11568; +constexpr double MassLambda0Bar = 1.11568; +constexpr double MassLambda1520 = 1.5195; +constexpr double MassSigmaMinus = 1.19744; +constexpr double MassSigmaBarPlus = 1.19744; +constexpr double MassSigmaPlus = 1.18937; +constexpr double MassSigmaBarMinus = 1.18937; +constexpr double MassSigma0 = 1.192642; +constexpr double MassSigma0Bar = 1.192642; +constexpr double MassXiMinus = 1.32171; +constexpr double MassXiPlusBar = 1.32171; +constexpr double MassOmegaMinus = 1.67245; +constexpr double MassOmegaPlusBar = 1.67245; + +// END OF THE GENERATED BLOCK + +// legacy names +constexpr double MassPhoton = MassGamma; +constexpr double MassMuon = MassMuonMinus; +constexpr double MassPionCharged = MassPiPlus; +constexpr double MassPionNeutral = MassPi0; +constexpr double MassKaonCharged = MassKPlus; +constexpr double MassKaonNeutral = MassK0; +constexpr double MassLambda = MassLambda0; +constexpr double MassHyperhydrog4 = MassHyperHydrogen4; +constexpr double MassHyperhelium4 = MassHyperHelium4; constexpr float LightSpeedCm2S = 299792458.e2; // C in cm/s constexpr float LightSpeedCm2NS = LightSpeedCm2S * 1e-9; // C in cm/ns -} // namespace physics -} // namespace constants -} // namespace o2 +} // namespace o2::constants::physics #endif diff --git a/Common/Constants/include/CommonConstants/make_pdg_header.py b/Common/Constants/include/CommonConstants/make_pdg_header.py new file mode 100755 index 0000000000000..409c6fe28b1c0 --- /dev/null +++ b/Common/Constants/include/CommonConstants/make_pdg_header.py @@ -0,0 +1,179 @@ +#!/usr/bin/env python3 + +# 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. + +"""! +@brief Generates the body of a C++ header with PDG codes and particle masses. +@author Vít Kučera , Inha University +@date 2023-09-21 +""" + +import os +from ctypes import c_bool +from enum import Enum + +import ROOT # pylint: disable=import-error + +name_script = os.path.basename(__file__) + +# Enum of PDG_t particles +class PdgROOT(Enum): + kDown = ROOT.kDown + kDownBar = ROOT.kDownBar + kUp = ROOT.kUp + kUpBar = ROOT.kUpBar + kStrange = ROOT.kStrange + kStrangeBar = ROOT.kStrangeBar + kCharm = ROOT.kCharm + kCharmBar = ROOT.kCharmBar + kBottom = ROOT.kBottom + kBottomBar = ROOT.kBottomBar + kTop = ROOT.kTop + kTopBar = ROOT.kTopBar + kGluon = ROOT.kGluon + kElectron = ROOT.kElectron + kPositron = ROOT.kPositron + kNuE = ROOT.kNuE + kNuEBar = ROOT.kNuEBar + kMuonMinus = ROOT.kMuonMinus + kMuonPlus = ROOT.kMuonPlus + kNuMu = ROOT.kNuMu + kNuMuBar = ROOT.kNuMuBar + kTauMinus = ROOT.kTauMinus + kTauPlus = ROOT.kTauPlus + kNuTau = ROOT.kNuTau + kNuTauBar = ROOT.kNuTauBar + kGamma = ROOT.kGamma + kZ0 = ROOT.kZ0 + kWPlus = ROOT.kWPlus + kWMinus = ROOT.kWMinus + kPi0 = ROOT.kPi0 + kK0Long = ROOT.kK0Long + kPiPlus = ROOT.kPiPlus + kPiMinus = ROOT.kPiMinus + kProton = ROOT.kProton + kProtonBar = ROOT.kProtonBar + kNeutron = ROOT.kNeutron + kNeutronBar = ROOT.kNeutronBar + kK0Short = ROOT.kK0Short + kK0 = ROOT.kK0 + kK0Bar = ROOT.kK0Bar + kKPlus = ROOT.kKPlus + kKMinus = ROOT.kKMinus + kLambda0 = ROOT.kLambda0 + kLambda0Bar = ROOT.kLambda0Bar + kLambda1520 = ROOT.kLambda1520 + kSigmaMinus = ROOT.kSigmaMinus + kSigmaBarPlus = ROOT.kSigmaBarPlus + kSigmaPlus = ROOT.kSigmaPlus + kSigmaBarMinus = ROOT.kSigmaBarMinus + kSigma0 = ROOT.kSigma0 + kSigma0Bar = ROOT.kSigma0Bar + kXiMinus = ROOT.kXiMinus + kXiPlusBar = ROOT.kXiPlusBar + kOmegaMinus = ROOT.kOmegaMinus + kOmegaPlusBar = ROOT.kOmegaPlusBar + + +# Enum of additional particles +class Pdg(Enum): + kB0 = 511 + kB0Bar = -511 + kBPlus = 521 + kBS = 531 + kBSBar = -531 + kD0 = 421 + kD0Bar = -421 + kDMinus = -411 + kDPlus = 411 + kDS = 431 + kDSBar = -431 + kDStar = 413 + kChiC1 = 20443 + kJPsi = 443 + kLambdaB0 = 5122 + kLambdaCPlus = 4122 + kOmegaC0 = 4332 + kPhi = 333 + kSigmaC0 = 4112 + kSigmaCPlusPlus = 4222 + kX3872 = 9920443 + kXi0 = 3322 + kXiB0 = 5232 + kXiCCPlusPlus = 4422 + kXiCPlus = 4232 + kXiCZero = 4132 + kDeuteron = 1000010020 + kTriton = 1000010030 + kHelium3 = 1000020030 + kAlpha = 1000020040 + kHyperTriton = 1010010030 + kHyperHydrogen4 = 1010010040 + kHyperHelium4 = 1010020040 + + +dbPdg = ROOT.o2.O2DatabasePDG + + +def mass(code): + """Returns particle mass from o2::O2DatabasePDG.""" + # Missing particles should be added in O2DatabasePDG.h. + success = c_bool(True) + return dbPdg.Mass(code, success) + + +def declare_mass(pdg, type="double") -> str: + """Returns a C++ declaration of a particle mass constant.""" + return f"constexpr {type} Mass{pdg.name[1:]} = {mass(pdg.value)};\n" + + +# Comment at the beginning of the output +str_block_begin = f"""// BEGINNING OF THE GENERATED BLOCK. +// DO NOT EDIT THIS BLOCK DIRECTLY! +// It has been generated by the {name_script} script. +// For modifications, edit the script and generate this block again. +""" +# Comment at the end of the output +str_block_end = """// END OF THE GENERATED BLOCK +""" +# Start of enum declarations of additional particles +str_enum_head = """/// \\brief Declarations of named PDG codes of particles missing in ROOT PDG_t +/// \\note Follow kCamelCase naming convention +/// \\link https://root.cern/doc/master/TPDGCode_8h.html +enum Pdg { +""" +# End of enum declarations of additional particles +str_enum_foot = "};\n" +# Documentation string for mass declarations of additional particles +str_mass_o2_head = """/// \\brief Declarations of masses for additional particles +""" +# Documentation string for mass declarations of PDG_t particles +str_mass_root_head = """/// \\brief Declarations of masses for particles in ROOT PDG_t +""" + +# Additional particles +str_enum = str_enum_head +str_mass_o2 = str_mass_o2_head +for c in Pdg: + str_enum += f" {c.name} = {c.value},\n" + str_mass_o2 += declare_mass(c) +str_enum = str_enum[:-2] + "\n" # Remove the last comma. +str_enum += str_enum_foot + +# PDG_t particles +str_mass_root = str_mass_root_head +for d in PdgROOT: + str_mass_root += declare_mass(d) + +# Header body +str_header = "\n".join((str_block_begin, str_enum, str_mass_o2, str_mass_root, str_block_end)) +print(str_header) diff --git a/Common/Field/src/MagFieldFast.cxx b/Common/Field/src/MagFieldFast.cxx index e652dffa396c5..5caad34d56dd4 100644 --- a/Common/Field/src/MagFieldFast.cxx +++ b/Common/Field/src/MagFieldFast.cxx @@ -238,7 +238,7 @@ bool MagFieldFast::Field(const math_utils::Point3D xyz, double bxyz[3]) bool MagFieldFast::GetSegment(float x, float y, float z, int& zSeg, int& rSeg, int& quadrant) const { // get segment of point location - const float zGridSpaceInv = 1.f / (kSolZMax * 2 / kNSolZRanges); + const float zGridSpaceInv = 1.f / (kSolZMax * 2 / (float)kNSolZRanges); zSeg = -1; if (z < kSolZMax) { if (z > -kSolZMax) { diff --git a/Common/Utils/include/CommonUtils/ConfigurableParam.h b/Common/Utils/include/CommonUtils/ConfigurableParam.h index 08356dc462de8..7099e37d5bc50 100644 --- a/Common/Utils/include/CommonUtils/ConfigurableParam.h +++ b/Common/Utils/include/CommonUtils/ConfigurableParam.h @@ -21,6 +21,7 @@ #include #include #include +#include class TFile; class TRootIOCtor; diff --git a/Common/Utils/include/CommonUtils/ConfigurationMacroHelper.h b/Common/Utils/include/CommonUtils/ConfigurationMacroHelper.h index be9bbcad25772..f16e501d96d40 100644 --- a/Common/Utils/include/CommonUtils/ConfigurationMacroHelper.h +++ b/Common/Utils/include/CommonUtils/ConfigurationMacroHelper.h @@ -69,21 +69,6 @@ T GetFromMacro(const std::string& file, const std::string& funcname, const std:: return *ptr; } -// just-in-time interpret some C++ function using ROOT and make result available to runtime -// functiondecl: A string coding the function to call (example "bool foo(){ return true; }") -// funcname: The name of the function to call (example "foo()") -// type: Return type of function (example "bool") -// unique: Some unique string identifier under which the result will be stored in gROOT global variable space -template -T JITAndEvalFunction(const std::string& functiondecl, const std::string& funcname, const std::string& type, const std::string& unique) -{ - /** interpret and execute a function and retrieve pointer to the returned type **/ - auto line = Form("%s; %s __%s__ = %s;", functiondecl.c_str(), type.c_str(), unique.c_str(), funcname.c_str()); - gROOT->ProcessLine(line); - auto ptr = (T*)gROOT->GetGlobal(Form("__%s__", unique.c_str()))->GetAddress(); - return *ptr; -} - } // namespace conf } // namespace o2 diff --git a/DataFormats/Detectors/HMPID/src/Cluster.cxx b/DataFormats/Detectors/HMPID/src/Cluster.cxx index 203609ddf8e1b..e14075c16c842 100644 --- a/DataFormats/Detectors/HMPID/src/Cluster.cxx +++ b/DataFormats/Detectors/HMPID/src/Cluster.cxx @@ -288,7 +288,7 @@ int Cluster::solve(std::vector* pCluLst, float* pSigmaCut, b } // Phase 0. Initialise Fitter - double arglist[10]; + double arglist[10]{0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; float ierflg = 0.; TVirtualFitter* fitter = TVirtualFitter::Fitter((TObject*)this, 3 * 6); // initialize Fitter if (fitter == nullptr) { @@ -349,10 +349,11 @@ int Cluster::solve(std::vector* pCluLst, float* pSigmaCut, b mSt = kMax; pCluLst->push_back(o2::hmpid::Cluster(*this)); //...add this raw cluster pCluLst->back().cleanPointers(); - } else { // or resonable number of local maxima to fit and user requested it + } else { // or resonable number of local maxima to fit and user requested it // Now ready for minimization step - arglist[0] = 500; // number of steps and sigma on pads charges - arglist[1] = 1.; // + arglist[0] = 500; // number of steps and sigma on pads charges + arglist[1] = 1.; // + /* ierflg = fitter->ExecuteCommand("SIMPLEX", arglist, 2); // start fitting with Simplex if (!ierflg) { fitter->ExecuteCommand("MIGRAD", arglist, 2); // fitting improved by Migrad @@ -367,6 +368,14 @@ int Cluster::solve(std::vector* pCluLst, float* pSigmaCut, b } } } + */ + + double strategy = 2.; + ierflg = fitter->ExecuteCommand("SET STR", &strategy, 1); // change level of strategy + if (!ierflg) { + fitter->ExecuteCommand("MIGRAD", arglist, 2); // fitting improved by Migrad + } + if (ierflg) { mSt = kAbn; // no convergence of the fit... } @@ -381,6 +390,7 @@ int Cluster::solve(std::vector* pCluLst, float* pSigmaCut, b fitter->GetParameter(3 * i + 1, sName, mYY, mErrY, dummy, dummy); // Y fitter->GetParameter(3 * i + 2, sName, mQ, mErrQ, dummy, dummy); // Q fitter->GetStats(mChi2, edm, errdef, nvpar, nparx); // get fit infos + // Printf("********************loc. max. = %i, X= %f, Y = %f, Q = %f**************************",i,mXX,mYY,mQ); if (mNlocMax > 1) { findClusterSize(i, pSigmaCut); // find clustersize for deconvoluted clusters // after this call, fSi temporarly is the calculated size. Later is set again @@ -401,6 +411,7 @@ int Cluster::solve(std::vector* pCluLst, float* pSigmaCut, b } } // setClusterParams(mXX, mYY, mCh); //need to fill the AliCluster3D part + // Printf("********************loc. max. = %i, X= %f, Y = %f, Q = %f**************************",i,mXX,mYY,mQ); pCluLst->push_back(o2::hmpid::Cluster(*this)); // add new unfolded cluster pCluLst->back().cleanPointers(); if (mNlocMax > 1) { @@ -438,12 +449,13 @@ Bool_t Cluster::isInPc() // Check if (X,Y) position is inside the PC limits // Arguments: // Returns: True or False + const auto param = o2::hmpid::Param::instanceNoGeo(); if (!mDigs) { LOGP(fatal, "digits are missing in the cluster"); } int pc = (*mDigs)[0]->getPh(); // (o2::hmpid::Digit*)&mDigs.at(iDig) - if (mXX < Param::minPcX(pc) || mXX > Param::maxPcX(pc) || mYY < Param::minPcY(pc) || mYY > Param::maxPcY(pc)) { + if (mXX < param->minPcX(pc) || mXX > param->maxPcX(pc) || mYY < param->minPcY(pc) || mYY > param->maxPcY(pc)) { return false; } diff --git a/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackParametrization.h b/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackParametrization.h index 880c2a59ac87f..105db7fbb4778 100644 --- a/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackParametrization.h +++ b/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackParametrization.h @@ -115,6 +115,8 @@ GPUconstexpr() int CovarMap[kNParams][kNParams] = {{0, 1, 3, 6, 10}, GPUconstexpr() int DiagMap[kNParams] = {0, 2, 5, 9, 14}; constexpr float HugeF = o2::constants::math::VeryBig; +constexpr float MaxPT = 100000.; // do not allow pTs exceeding this value (to avoid NANs) +constexpr float MinPTInv = 1. / MaxPT; // do not allow q/pTs less this value (to avoid NANs) template class TrackParametrization @@ -541,7 +543,10 @@ template GPUdi() auto TrackParametrization::getPtInv() const -> value_t { // return the inverted track pT - const value_t ptInv = gpu::CAMath::Abs(mP[kQ2Pt]); + value_t ptInv = gpu::CAMath::Abs(mP[kQ2Pt]); + if (ptInv < MinPTInv) { + ptInv = MinPTInv; + } return (mAbsCharge > 1) ? ptInv / mAbsCharge : ptInv; } @@ -550,8 +555,8 @@ template GPUdi() auto TrackParametrization::getP2Inv() const -> value_t { // return the inverted track momentum^2 - const value_t p2 = mP[kQ2Pt] * mP[kQ2Pt] / (1.f + getTgl() * getTgl()); - return (mAbsCharge > 1) ? p2 / (mAbsCharge * mAbsCharge) : p2; + value_t p2 = getPtInv(); + return p2 * p2 / (1.f + getTgl() * getTgl()); } //____________________________________________________________ @@ -559,17 +564,15 @@ template GPUdi() auto TrackParametrization::getP2() const -> value_t { // return the track momentum^2 - const value_t p2inv = getP2Inv(); - return (p2inv > o2::constants::math::Almost0) ? 1.f / p2inv : o2::constants::math::VeryBig; + return 1.f / getP2Inv(); // getP2Inv is protected against being 0, full charge accounted } //____________________________________________________________ template GPUdi() auto TrackParametrization::getPInv() const -> value_t { - // return the inverted track momentum^2 - const value_t pInv = gpu::CAMath::Abs(mP[kQ2Pt]) / gpu::CAMath::Sqrt(1.f + getTgl() * getTgl()); - return (mAbsCharge > 1) ? pInv / mAbsCharge : pInv; + // return the inverted track momentum + return getPtInv() / gpu::CAMath::Sqrt(1.f + getTgl() * getTgl()); // getPtInv() is protected against being 0, full charge accounted } //____________________________________________________________ @@ -577,8 +580,7 @@ template GPUdi() auto TrackParametrization::getP() const -> value_t { // return the track momentum - const value_t pInv = getPInv(); - return (pInv > o2::constants::math::Almost0) ? 1.f / pInv : o2::constants::math::VeryBig; + return 1.f / getPInv(); // getPInv is already protected against being 0 } //____________________________________________________________ @@ -586,11 +588,7 @@ template GPUdi() auto TrackParametrization::getPt() const -> value_t { // return the track transverse momentum - value_t ptI = gpu::CAMath::Abs(mP[kQ2Pt]); - if (mAbsCharge > 1) { - ptI /= mAbsCharge; - } - return (ptI > o2::constants::math::Almost0) ? 1.f / ptI : o2::constants::math::VeryBig; + return 1.f / getPtInv(); // getPtInv is already protected against being 0 } //____________________________________________________________ diff --git a/DataFormats/Reconstruction/src/TrackParametrization.cxx b/DataFormats/Reconstruction/src/TrackParametrization.cxx index e084315a0d0b1..b0f8febf3ba59 100644 --- a/DataFormats/Reconstruction/src/TrackParametrization.cxx +++ b/DataFormats/Reconstruction/src/TrackParametrization.cxx @@ -144,9 +144,9 @@ template GPUd() bool TrackParametrization::getPosDirGlo(gpu::gpustd::array& posdirp) const { // fill vector with lab x,y,z,px/p,py/p,pz/p,p,sinAlpha,cosAlpha - value_t ptI = gpu::CAMath::Abs(getQ2Pt()); + value_t ptI = getPtInv(); value_t snp = getSnp(); - if (ptI < constants::math::Almost0 || gpu::CAMath::Abs(snp) > constants::math::Almost1) { + if (gpu::CAMath::Abs(snp) > constants::math::Almost1) { return false; } value_t &sn = posdirp[7], &cs = posdirp[8]; diff --git a/Detectors/Base/CMakeLists.txt b/Detectors/Base/CMakeLists.txt index c7ea137679fb1..8c1544164a851 100644 --- a/Detectors/Base/CMakeLists.txt +++ b/Detectors/Base/CMakeLists.txt @@ -25,7 +25,6 @@ o2_add_library(DetectorsBase src/TFIDInfoHelper.cxx src/GRPGeomHelper.cxx src/Stack.cxx - src/VMCSeederService.cxx PUBLIC_LINK_LIBRARIES FairRoot::Base O2::CommonUtils O2::DetectorsCommonDataFormats diff --git a/Detectors/Base/include/DetectorsBase/Propagator.h b/Detectors/Base/include/DetectorsBase/Propagator.h index 2227f2156c401..600e71e7fb95d 100644 --- a/Detectors/Base/include/DetectorsBase/Propagator.h +++ b/Detectors/Base/include/DetectorsBase/Propagator.h @@ -137,6 +137,7 @@ class PropagatorImpl GPUd() bool hasMagFieldSet() const { return mField != nullptr; } GPUd() void estimateLTFast(o2::track::TrackLTIntegral& lt, const o2::track::TrackParametrization& trc) const; + GPUd() float estimateLTIncrement(const o2::track::TrackParametrization& trc, const o2::math_utils::Point3D& postStart, const o2::math_utils::Point3D& posEnd) const; #ifndef GPUCA_GPUCODE static PropagatorImpl* Instance(bool uninitialized = false) diff --git a/Detectors/Base/include/DetectorsBase/VMCSeederService.h b/Detectors/Base/include/DetectorsBase/VMCSeederService.h deleted file mode 100644 index 1669c73b39620..0000000000000 --- a/Detectors/Base/include/DetectorsBase/VMCSeederService.h +++ /dev/null @@ -1,50 +0,0 @@ -// 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 ALICEO2_VMC_SEEDERSERVICE_H -#define ALICEO2_VMC_SEEDERSERVICE_H - -class TVirtualMC; - -#include - -namespace o2 -{ -namespace base -{ - -// A simple (singleton) service to propagate random number seeding -// to VMC engines in a loosely-coupled, just-in-time way (without needing to compile against these libraries). -class VMCSeederService -{ - - public: - static VMCSeederService const& instance() - { - static VMCSeederService inst; - return inst; - } - - void setSeed() const; // will propagate seed to the VMC engines - - typedef std::function SeederFcn; - - private: - VMCSeederService(); - void initSeederFunction(TVirtualMC const*); - - SeederFcn mSeederFcn; // the just-in-time compiled function talking to the VMC engines -}; - -} // namespace base -} // namespace o2 - -#endif diff --git a/Detectors/Base/src/Propagator.cxx b/Detectors/Base/src/Propagator.cxx index cdc301a47cb6b..54baabb8f3af5 100644 --- a/Detectors/Base/src/Propagator.cxx +++ b/Detectors/Base/src/Propagator.cxx @@ -616,6 +616,27 @@ GPUd() bool PropagatorImpl::propagateToDCABxByBz(const math_utils::Poin return true; } +//____________________________________________________________ +template +GPUd() float PropagatorImpl::estimateLTIncrement(const o2::track::TrackParametrization& trc, + const o2::math_utils::Point3D& posStart, + const o2::math_utils::Point3D& posEnd) const +{ + // estimate helical step increment between 2 point + float dX = posEnd.X() - posStart.X(), dY = posEnd.Y() - posStart.Y(), dZ = posEnd.Z() - posStart.Z(), d2XY = dX * dX + dY * dY; + if (getNominalBz() != 0) { // circular arc = 2*R*asin(dXY/2R) + float b[3]; + o2::math_utils::Point3D posAv(0.5 * (posEnd.X() + posStart.X()), 0.5 * (posEnd.Y() + posStart.Y()), 0.5 * (posEnd.Z() + posStart.Z())); + getFieldXYZ(posAv, b); + float curvH = math_utils::detail::abs(0.5f * trc.getCurvature(b[2])), asArg = curvH * math_utils::detail::sqrt(d2XY); + if (curvH > 0.f) { + d2XY = asArg < 1.f ? math_utils::detail::asin(asArg) / curvH : o2::constants::math::PIHalf / curvH; + d2XY *= d2XY; + } + } + return math_utils::detail::sqrt(d2XY + dZ * dZ); +} + //____________________________________________________________ template GPUd() void PropagatorImpl::estimateLTFast(o2::track::TrackLTIntegral& lt, const o2::track::TrackParametrization& trc) const diff --git a/Detectors/Base/src/Stack.cxx b/Detectors/Base/src/Stack.cxx index 64daf2eccefc5..2e3de9b75bffc 100644 --- a/Detectors/Base/src/Stack.cxx +++ b/Detectors/Base/src/Stack.cxx @@ -15,7 +15,6 @@ #include "DetectorsBase/Stack.h" #include "DetectorsBase/Detector.h" -#include #include "DetectorsCommonDataFormats/DetID.h" #include "SimulationDataFormat/MCTrack.h" @@ -363,7 +362,6 @@ TParticle* Stack::PopNextTrack(Int_t& iTrack) // LOG(info) << "SEEDING NEW TRACK USING HASH" << hash; // init seed per track gRandom->SetSeed(hash); - o2::base::VMCSeederService::instance().setSeed(); // NOTE: THE BETTER PLACE WOULD BE IN PRETRACK HOOK BUT THIS DOES NOT SEEM TO WORK // WORKS ONLY WITH G3 SINCE G4 DOES NOT CALL THIS FUNCTION } // .trackSeed ? diff --git a/Detectors/Base/src/VMCSeederService.cxx b/Detectors/Base/src/VMCSeederService.cxx deleted file mode 100644 index 95dc90da6efdf..0000000000000 --- a/Detectors/Base/src/VMCSeederService.cxx +++ /dev/null @@ -1,49 +0,0 @@ -// Copyright 2019-2020 CERN and copyright holders of ALICE O2. -// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. -// All rights not expressly granted are reserved. -// -// This software is distributed under the terms of the GNU General Public -// License v3 (GPL Version 3), copied verbatim in the file "COPYING". -// -// In applying this license CERN does not waive the privileges and immunities -// granted to it by virtue of its status as an Intergovernmental Organization -// or submit itself to any jurisdiction. - -#include "DetectorsBase/VMCSeederService.h" -#include "TVirtualMC.h" -#include // for FairLogger -#include // for ROOT JIT helpers - -using namespace o2::base; - -void VMCSeederService::initSeederFunction(TVirtualMC const* vmc) -{ - if (strcmp(vmc->GetName(), "TGeant3TGeo") == 0) { - // Geant3 doesn't need anything special in our context - mSeederFcn = []() {}; - } else if (strcmp(vmc->GetName(), "TGeant4") == 0) { - // dynamically get access to the Geant4_VMC seeding function (without this function linking against Geant4) - std::string G4func("std::function G4func() { gSystem->Load(\"libgeant4vmc\"); return [](){ ((TGeant4*)TVirtualMC::GetMC())->SetRandomSeed(); };}"); - mSeederFcn = o2::conf::JITAndEvalFunction(G4func, "G4func()", "std::function", "VMCSEEDERFUNC123"); - } else { - LOG(warn) << "Unknown VMC engine or unimplemented VMC seeding function"; - mSeederFcn = []() {}; - } -} - -// constructor -VMCSeederService::VMCSeederService() -{ - auto vmc = TVirtualMC::GetMC(); - if (vmc) { - LOG(info) << "Seeder initializing for " << vmc->GetName(); - initSeederFunction(vmc); - } else { - LOG(fatal) << " Seeder could not be initialized (no VMC instance found)"; - } -} - -void VMCSeederService::setSeed() const -{ - mSeederFcn(); -} diff --git a/Detectors/FIT/workflow/include/FITWorkflow/FITIntegrateClusterSpec.h b/Detectors/FIT/workflow/include/FITWorkflow/FITIntegrateClusterSpec.h index c46c07a4238a7..603c123d337c0 100644 --- a/Detectors/FIT/workflow/include/FITWorkflow/FITIntegrateClusterSpec.h +++ b/Detectors/FIT/workflow/include/FITWorkflow/FITIntegrateClusterSpec.h @@ -180,9 +180,9 @@ o2::framework::DataProcessorSpec getFITIntegrateClusterSpec(const bool disableWr inputs); std::vector outputs; - outputs.emplace_back(FitType::getDataOrigin(), FitType::getDataDescriptionFITC(), 0, Lifetime::Sporadic); + outputs.emplace_back(FitType::getDataOrigin(), FitType::getDataDescriptionFITC(), 0, Lifetime::Timeframe); if (!disableWriter) { - outputs.emplace_back(FitType::getDataOrigin(), FitType::getDataDescriptionFITTFId(), 0, Lifetime::Sporadic); + outputs.emplace_back(FitType::getDataOrigin(), FitType::getDataDescriptionFITTFId(), 0, Lifetime::Timeframe); } return DataProcessorSpec{ diff --git a/Detectors/FIT/workflow/include/FITWorkflow/FITIntegrateClusterWriterSpec.h b/Detectors/FIT/workflow/include/FITWorkflow/FITIntegrateClusterWriterSpec.h index c5d0afe798740..c53d3e9b6425c 100644 --- a/Detectors/FIT/workflow/include/FITWorkflow/FITIntegrateClusterWriterSpec.h +++ b/Detectors/FIT/workflow/include/FITWorkflow/FITIntegrateClusterWriterSpec.h @@ -38,8 +38,8 @@ DataProcessorSpec getFITIntegrateClusterWriterSpec() return MakeRootTreeWriterSpec(fmt::format("{}-currents-writer", FitType::getName()).data(), treeFile.data(), treeName.data(), - BranchDefinition::DataTStruct>{InputSpec{"ifitc", FitType::getDataOrigin(), FitType::getDataDescriptionFITC(), 0}, "IFITC", 1}, - BranchDefinition{InputSpec{"tfID", FitType::getDataOrigin(), FitType::getDataDescriptionFITTFId(), 0}, "tfID", 1})(); + BranchDefinition::DataTStruct>{InputSpec{"ifitc", FitType::getDataOrigin(), FitType::getDataDescriptionFITC(), 0, Lifetime::Timeframe}, "IFITC", 1}, + BranchDefinition{InputSpec{"tfID", FitType::getDataOrigin(), FitType::getDataDescriptionFITTFId(), 0, Lifetime::Timeframe}, "tfID", 1})(); } } // end namespace fit diff --git a/Detectors/GlobalTracking/include/GlobalTracking/MatchTOF.h b/Detectors/GlobalTracking/include/GlobalTracking/MatchTOF.h index 607ad27eb3c20..d831a9c71d1fb 100644 --- a/Detectors/GlobalTracking/include/GlobalTracking/MatchTOF.h +++ b/Detectors/GlobalTracking/include/GlobalTracking/MatchTOF.h @@ -154,24 +154,6 @@ class MatchTOF void setTPCVDrift(const o2::tpc::VDriftCorrFact& v); void setTPCCorrMaps(o2::gpu::CorrectionMapsHelper* maph); - ///< set input TPC tracks cluster indices - void setTPCTrackClusIdxInp(const gsl::span inp) - { - mTPCTrackClusIdx = inp; - } - - ///< set input TPC cluster sharing map - void setTPCClustersSharingMap(const gsl::span inp) - { - mTPCRefitterShMap = inp; - } - - ///< set input TPC clusters - void setTPCClustersInp(const o2::tpc::ClusterNativeAccess* inp) - { - mTPCClusterIdxStruct = inp; - } - void setFIT(bool value = true) { mIsFIT = value; } static int findFITIndex(int bc, const gsl::span& FITRecPoints, unsigned long firstOrbit); @@ -179,18 +161,32 @@ class MatchTOF bool makeConstrainedTPCTrack(int matchedID, o2::dataformats::TrackTPCTOF& trConstr); ///< populate externally provided container by TOF-time-constrained TPC tracks - template - void makeConstrainedTPCTracks(V& container) + template + void makeConstrainedTPCTracks(MtcInfo& mtcCont, MCInfo& MCCont, CTrack& trcCont) { + auto nmatch = getMatchedTrackVector(o2::dataformats::MatchInfoTOFReco::TrackType::TPC).size(); // preliminary matches + mtcCont.reserve(nmatch); + trcCont.reserve(nmatch); + if (mMCTruthON) { + MCCont.reserve(nmatch); + } checkRefitter(); - int nmatched = mMatchedTracks[trkType::TPC].size(), nconstrained = 0; - container.resize(nmatched); - for (unsigned i = 0; i < nmatched; i++) { - if (makeConstrainedTPCTrack(i, container[nconstrained])) { + + auto& mclabs = getMatchedTOFLabelsVector(o2::dataformats::MatchInfoTOFReco::TrackType::TPC); + auto& info = getMatchedTrackVector(o2::dataformats::MatchInfoTOFReco::TrackType::TPC); + int nconstrained = 0; + for (unsigned i = 0; i < nmatch; i++) { + auto& ctr = trcCont.emplace_back(); + if (makeConstrainedTPCTrack(i, ctr)) { + mtcCont.push_back(info[i]); + if (mMCTruthON) { + MCCont.push_back(mclabs[i]); + } nconstrained++; + } else { + trcCont.pop_back(); } } - container.resize(nconstrained); } void setTS(unsigned long creationTime) { mTimestamp = creationTime; } diff --git a/Detectors/GlobalTracking/src/MatchHMP.cxx b/Detectors/GlobalTracking/src/MatchHMP.cxx index a6706d6408624..207bbdb9bfbf6 100644 --- a/Detectors/GlobalTracking/src/MatchHMP.cxx +++ b/Detectors/GlobalTracking/src/MatchHMP.cxx @@ -524,7 +524,7 @@ void MatchHMP::doMatching() // 6. Set match information int cluSize = bestHmpCluster->size(); - matching.setHMPIDmip(bestHmpCluster->x(), bestHmpCluster->y(), bestHmpCluster->q(), 0); // store mip info in any case + matching.setHMPIDmip(bestHmpCluster->x(), bestHmpCluster->y(), (int)bestHmpCluster->q(), 0); // store mip info in any case matching.setMipClusSize(bestHmpCluster->size()); matching.setIdxHMPClus(iCh, index + 1000 * cluSize); // set chamber, index of cluster + cluster size matching.setHMPIDtrk(xPc, yPc, theta, phi); diff --git a/Detectors/GlobalTracking/src/MatchTOF.cxx b/Detectors/GlobalTracking/src/MatchTOF.cxx index 0eb9cf430beec..4b38f59814422 100644 --- a/Detectors/GlobalTracking/src/MatchTOF.cxx +++ b/Detectors/GlobalTracking/src/MatchTOF.cxx @@ -484,6 +484,13 @@ bool MatchTOF::prepareTPCData() } // loop over tracks of single sector } + if (mRecoCont->inputsTPCclusters) { + mTPCClusterIdxStruct = &mRecoCont->inputsTPCclusters->clusterIndex; + mTPCTracksArray = mRecoCont->getTPCTracks(); + mTPCTrackClusIdx = mRecoCont->getTPCTracksClusterRefs(); + mTPCRefitterShMap = mRecoCont->clusterShMapTPC; + } + return true; } @@ -1752,10 +1759,10 @@ bool MatchTOF::makeConstrainedTPCTrack(int matchedID, o2::dataformats::TrackTPCT float chi2 = 0; mTPCRefitter->setTrackReferenceX(o2::constants::geom::XTPCInnerRef); if (mTPCRefitter->RefitTrackAsTrackParCov(trConstr, tpcTrOrig.getClusterRef(), timeTOFTB, &chi2, false, true) < 0) { // outward refit after resetting cov.mat. - LOG(debug) << "Inward Refit failed"; + LOGP(debug, "Inward Refit failed {}", trConstr.asString()); return false; } - const auto& trackTune = o2::globaltracking::TrackTuneParams::Instance(); // if needed, correct TPC track after inwar refit + const auto& trackTune = o2::globaltracking::TrackTuneParams::Instance(); // if needed, correct TPC track after inward refit if (!trackTune.useTPCInnerCorr) { trConstr.updateParams(trackTune.tpcParInner); } @@ -1767,7 +1774,7 @@ bool MatchTOF::makeConstrainedTPCTrack(int matchedID, o2::dataformats::TrackTPCT // mTPCRefitter->setTrackReferenceX(o2::constants::geom::XTPCOuterRef); if (mTPCRefitter->RefitTrackAsTrackParCov(trConstrOut, tpcTrOrig.getClusterRef(), timeTOFTB, &chi2, true, true) < 0) { // outward refit after resetting cov.mat. - LOG(debug) << "Outward refit failed"; + LOGP(debug, "Outward Refit failed {}", trConstrOut.asString()); return false; } } diff --git a/Detectors/GlobalTracking/src/MatchTPCITS.cxx b/Detectors/GlobalTracking/src/MatchTPCITS.cxx index dda53b2598d6e..8d5350c2da106 100644 --- a/Detectors/GlobalTracking/src/MatchTPCITS.cxx +++ b/Detectors/GlobalTracking/src/MatchTPCITS.cxx @@ -1442,16 +1442,7 @@ bool MatchTPCITS::refitTrackTPCITS(int iTPC, int& iITS, pmr::vector posAv(0.5 * (posEnd.x() + posStart.x()), 0.5 * (posEnd.y() + posStart.y()), 0.5 * (posEnd.z() + posStart.z())); - propagator->getFieldXYZ(posAv, b); - float curvH = std::abs(0.5f * tracOut.getCurvature(b[2])), arcXY = 1. / curvH * std::asin(curvH * std::sqrt(d2XY)); - d2XY = arcXY * arcXY; - } - auto lInt = std::sqrt(d2XY + dZ * dZ); + auto lInt = propagator->estimateLTIncrement(tracOut, posStart, posEnd); tofL.addStep(lInt, tracOut.getP2Inv()); tofL.addX2X0(lInt * mTPCmeanX0Inv); propagator->PropagateToXBxByBz(tracOut, o2::constants::geom::XTPCOuterRef, MaxSnp, 10., mUseMatCorrFlag, &tofL); @@ -1464,7 +1455,6 @@ bool MatchTPCITS::refitTrackTPCITS(int iTPC, int& iITS, pmr::vectorestimateLTFast(tofL, winLink); // guess about initial value for the track integral from the origin - // refit track outward in the ITS const auto& itsClRefs = ABTrackletRefs[iITSAB]; int nclRefit = 0, ncl = itsClRefs.getNClusters(); @@ -1588,20 +1577,10 @@ bool MatchTPCITS::refitABTrack(int iITSAB, const TPCABSeed& seed, pmr::vector posAv(0.5 * (posEnd.x() + posStart.x()), 0.5 * (posEnd.y() + posStart.y()), 0.5 * (posEnd.z() + posStart.z())); - propagator->getFieldXYZ(posAv, b); - float curvH = std::abs(0.5f * tracOut.getCurvature(b[2])), arcXY = 1. / curvH * std::asin(curvH * std::sqrt(d2XY)); - d2XY = arcXY * arcXY; - } - auto lInt = std::sqrt(d2XY + dZ * dZ); + auto lInt = propagator->estimateLTIncrement(tracOut, posStart, posEnd); tofL.addStep(lInt, tracOut.getP2Inv()); tofL.addX2X0(lInt * mTPCmeanX0Inv); propagator->PropagateToXBxByBz(tracOut, o2::constants::geom::XTPCOuterRef, MaxSnp, 10., mUseMatCorrFlag, &tofL); - const auto& trackTune = o2::globaltracking::TrackTuneParams::Instance(); if (trackTune.useTPCOuterCorr) { tracOut.updateParams(trackTune.tpcParOuter); @@ -1888,6 +1867,7 @@ void MatchTPCITS::refitABWinners(pmr::vector& matc } lID = winL.parentID; } + clref.setEntries(ncl); if (!refitABTrack(ABTrackletRefs.size() - 1, ABSeed, matchedTracks, ABTrackletClusterIDs, ABTrackletRefs)) { // on failure, destroy added tracklet reference ABTrackletRefs.pop_back(); ABTrackletClusterIDs.resize(start); // RSS @@ -1896,7 +1876,6 @@ void MatchTPCITS::refitABWinners(pmr::vector& matc } continue; } - clref.setEntries(ncl); if (mMCTruthON) { o2::MCCompLabel lab; int maxL = 0; // find most encountered label diff --git a/Detectors/GlobalTrackingWorkflow/src/TOFMatcherSpec.cxx b/Detectors/GlobalTrackingWorkflow/src/TOFMatcherSpec.cxx index 0c8fd7a809a1c..134ffbb75f341 100644 --- a/Detectors/GlobalTrackingWorkflow/src/TOFMatcherSpec.cxx +++ b/Detectors/GlobalTrackingWorkflow/src/TOFMatcherSpec.cxx @@ -168,21 +168,15 @@ void TOFMatcherSpec::run(ProcessingContext& pc) mMatcher.setTS(creationTime); mMatcher.run(recoData, pc.services().get().firstTForbit); + static pmr::vector dummyMCLab; if (isTPCused) { - pc.outputs().snapshot(Output{o2::header::gDataOriginTOF, "MTC_TPC", ss, Lifetime::Timeframe}, mMatcher.getMatchedTrackVector(o2::dataformats::MatchInfoTOFReco::TrackType::TPC)); - if (mUseMC) { - pc.outputs().snapshot(Output{o2::header::gDataOriginTOF, "MCMTC_TPC", ss, Lifetime::Timeframe}, mMatcher.getMatchedTOFLabelsVector(o2::dataformats::MatchInfoTOFReco::TrackType::TPC)); - } - + auto& mtcInfo = pc.outputs().make>(Output{o2::header::gDataOriginTOF, "MTC_TPC", ss, Lifetime::Timeframe}); + auto& mclabels = mUseMC ? pc.outputs().make>(Output{o2::header::gDataOriginTOF, "MCMTC_TPC", ss, Lifetime::Timeframe}) : dummyMCLab; + auto& tracksTPCTOF = pc.outputs().make>(OutputRef{"tpctofTracks", ss}); auto nmatch = mMatcher.getMatchedTrackVector(o2::dataformats::MatchInfoTOFReco::TrackType::TPC).size(); - if (mDoTPCRefit) { - LOG(debug) << "Refitting " << nmatch << " matched TPC tracks with TOF time info"; - } else { - LOG(debug) << "Shifting Z for " << nmatch << " matched TPC tracks according to TOF time info"; - } - auto& tracksTPCTOF = pc.outputs().make>(OutputRef{"tpctofTracks", ss}, nmatch); - mMatcher.makeConstrainedTPCTracks(tracksTPCTOF); + LOG(debug) << (mDoTPCRefit ? "Refitting " : "Shifting Z for ") << nmatch << " matched TPC tracks with TOF time info"; + mMatcher.makeConstrainedTPCTracks(mtcInfo, mclabels, tracksTPCTOF); } if (isITSTPCused) { @@ -249,6 +243,9 @@ DataProcessorSpec getTOFMatcherSpec(GID::mask_t src, bool useMC, bool useFIT, bo } dataRequest->requestTracks(src, useMC); dataRequest->requestClusters(GID::getSourceMask(GID::TOF), useMC); + if (tpcRefit && src[GID::TPC]) { + dataRequest->requestClusters(GID::getSourceMask(GID::TPC), false); + } if (useFIT) { dataRequest->requestFT0RecPoints(false); } diff --git a/Detectors/GlobalTrackingWorkflow/src/tof-matcher-workflow.cxx b/Detectors/GlobalTrackingWorkflow/src/tof-matcher-workflow.cxx index 964e2d118c069..c63e35deda2ed 100644 --- a/Detectors/GlobalTrackingWorkflow/src/tof-matcher-workflow.cxx +++ b/Detectors/GlobalTrackingWorkflow/src/tof-matcher-workflow.cxx @@ -65,6 +65,7 @@ void customize(std::vector& workflowOptions) {"output-type", o2::framework::VariantType::String, "matching-info", {"matching-info, calib-info"}}, {"enable-dia", o2::framework::VariantType::Bool, false, {"to require diagnostic freq and then write to calib outputs (obsolete since now default)"}}, {"trd-extra-tolerance", o2::framework::VariantType::Float, 500.0f, {"Extra time tolerance for TRD tracks in ns"}}, + {"refit-tpc-tof", o2::framework::VariantType::Bool, false, {"Refit unconstrained TPC tracks matched to TOF (if false - just move)"}}, {"write-matchable", o2::framework::VariantType::Bool, false, {"write all matchable pairs in a file (o2matchable_tof.root)"}}, {"configKeyValues", VariantType::String, "", {"Semicolon separated key=value strings ..."}}, {"combine-devices", o2::framework::VariantType::Bool, false, {"merge DPL source/writer devices"}}}; @@ -98,6 +99,7 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) auto sclOpt = o2::tpc::CorrectionMapsLoader::parseGlobalOptions(configcontext.options()); bool writematching = 0; bool writecalib = 0; + bool refitTPCTOF = configcontext.options().get("refit-tpc-tof"); auto outputType = configcontext.options().get("output-type"); auto nLanes = configcontext.options().get("tof-lanes"); if (outputType.rfind("matching-info") < outputType.size()) { @@ -142,6 +144,9 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) if (useFIT) { clustermask |= GID::getSourceMask(GID::FT0); } + if (src[GID::TPC] && refitTPCTOF) { // load clusters + clustermask |= GID::getSourceMask(GID::TPC); + } if (sclOpt.lumiType == 1) { src = src | GID::getSourcesMask("CTP"); } @@ -163,7 +168,7 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) } } - specs.emplace_back(o2::globaltracking::getTOFMatcherSpec(src, useMC, useFIT, false, strict, extratolerancetrd, writeMatchable, sclOpt, nLanes)); // doTPCrefit not yet supported (need to load TPC clusters?) + specs.emplace_back(o2::globaltracking::getTOFMatcherSpec(src, useMC, useFIT, refitTPCTOF, strict, extratolerancetrd, writeMatchable, sclOpt, nLanes)); // doTPCrefit not yet supported (need to load TPC clusters?) if (!disableRootOut) { std::vector writers; diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/MathUtils.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/MathUtils.h index b9a2b599f0ba8..9093609144283 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/MathUtils.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/MathUtils.h @@ -77,14 +77,20 @@ GPUhdi() float math_utils::computeCurvature(float x1, float y1, float x2, float 0.5f * ((y3 - y2) * (y2 * y2 - y1 * y1 + x2 * x2 - x1 * x1) - (y2 - y1) * (y3 * y3 - y2 * y2 + x3 * x3 - x2 * x2)); const float b = 0.5f * ((x2 - x1) * (y3 * y3 - y2 * y2 + x3 * x3 - x2 * x2) - (x3 - x2) * (y2 * y2 - y1 * y1 + x2 * x2 - x1 * x1)); - - return -1.f * d / o2::gpu::CAMath::Sqrt((d * x1 - a) * (d * x1 - a) + (d * y1 - b) * (d * y1 - b)); + const float den2 = (d * x1 - a) * (d * x1 - a) + (d * y1 - b) * (d * y1 - b); + return den2 > 0.f ? -1.f * d / o2::gpu::CAMath::Sqrt(den2) : 0.f; } GPUhdi() float math_utils::computeCurvatureCentreX(float x1, float y1, float x2, float y2, float x3, float y3) { - const float k1 = (y2 - y1) / (x2 - x1), k2 = (y3 - y2) / (x3 - x2); - return 0.5f * (k1 * k2 * (y1 - y3) + k2 * (x1 + x2) - k1 * (x2 + x3)) / (k2 - k1); + float dx21 = x2 - x1, dx32 = x3 - x2; + if (dx21 == 0.f || dx32 == 0.f) { // add small offset + x2 += 1e-4; + dx21 = x2 - x1; + dx32 = x3 - x2; + } + float k1 = (y2 - y1) / dx21, k2 = (y3 - y2) / dx32; + return (k1 != k2) ? 0.5f * (k1 * k2 * (y1 - y3) + k2 * (x1 + x2) - k1 * (x2 + x3)) / (k2 - k1) : 1e5; } GPUhdi() float math_utils::computeTanDipAngle(float x1, float y1, float x2, float y2, float z1, float z2) diff --git a/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx b/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx index d5911c5a0f7e2..22e67cfe63461 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx @@ -961,16 +961,16 @@ track::TrackParCov TrackerTraits::buildTrackSeed(const Cluster& cluster1, const const float snp = zeroField ? tgp / std::sqrt(1.f + tgp * tgp) : crv * (x3 - math_utils::computeCurvatureCentreX(x3, y3, x2, y2, x1, y1)); const float tgl12 = math_utils::computeTanDipAngle(x1, y1, x2, y2, z1, z2); const float tgl23 = math_utils::computeTanDipAngle(x2, y2, x3, y3, z2, z3); - + const float q2pt = zeroField ? 1.f / o2::track::kMostProbablePt : crv / (getBz() * o2::constants::math::B2C); + const float q2pt2 = crv * crv; + const float sg2q2pt = track::kC1Pt2max * (q2pt2 > 0.0005 ? (q2pt2 < 1 ? q2pt2 : 1) : 0.0005); return track::TrackParCov(tf3.xTrackingFrame, tf3.alphaTrackingFrame, - {y3, z3, snp, 0.5f * (tgl12 + tgl23), - zeroField ? 1.f / o2::track::kMostProbablePt - : crv / (getBz() * o2::constants::math::B2C)}, + {y3, z3, snp, 0.5f * (tgl12 + tgl23), q2pt}, {tf3.covarianceTrackingFrame[0], tf3.covarianceTrackingFrame[1], tf3.covarianceTrackingFrame[2], 0.f, 0.f, track::kCSnp2max, 0.f, 0.f, 0.f, track::kCTgl2max, - 0.f, 0.f, 0.f, 0.f, track::kC1Pt2max}); + 0.f, 0.f, 0.f, 0.f, sg2q2pt}); } void TrackerTraits::setBz(float bz) diff --git a/Detectors/MUON/MCH/GlobalMapping/src/testChannelCode.cxx b/Detectors/MUON/MCH/GlobalMapping/src/testChannelCode.cxx index 8e6d777f8f492..68de2c4169f2e 100644 --- a/Detectors/MUON/MCH/GlobalMapping/src/testChannelCode.cxx +++ b/Detectors/MUON/MCH/GlobalMapping/src/testChannelCode.cxx @@ -17,6 +17,7 @@ #include #include +#include #include "MCHGlobalMapping/ChannelCode.h" BOOST_AUTO_TEST_CASE(CtorShowThrowForInvalidDeId) diff --git a/Detectors/MUON/MCH/Tracking/src/TrackFinderSpec.cxx b/Detectors/MUON/MCH/Tracking/src/TrackFinderSpec.cxx index c40df5d4b448e..6239186309dc3 100644 --- a/Detectors/MUON/MCH/Tracking/src/TrackFinderSpec.cxx +++ b/Detectors/MUON/MCH/Tracking/src/TrackFinderSpec.cxx @@ -66,7 +66,7 @@ class TrackFinderTask { public: //_________________________________________________________________________________________________ - TrackFinderTask(bool computeTime, bool digits, std::shared_ptr req) + TrackFinderTask(bool computeTime, bool digits, std::shared_ptr req) : mComputeTime(computeTime), mDigits(digits), mCCDBRequest(req) {} //_________________________________________________________________________________________________ diff --git a/Detectors/MUON/MID/Workflow/src/RawCheckerSpec.cxx b/Detectors/MUON/MID/Workflow/src/RawCheckerSpec.cxx index c2e9e3fc5c635..58774a6336294 100644 --- a/Detectors/MUON/MID/Workflow/src/RawCheckerSpec.cxx +++ b/Detectors/MUON/MID/Workflow/src/RawCheckerSpec.cxx @@ -48,7 +48,7 @@ template class RawCheckerDeviceDPL { public: - RawCheckerDeviceDPL(const std::vector& feeIds, const CrateMasks& crateMasks, const ElectronicsDelay& electronicsDelay) : mFeeIds(feeIds), mCrateMasks(crateMasks), mElectronicsDelay(electronicsDelay) {} + RawCheckerDeviceDPL(const std::vector& feeIds, const CrateMasks& crateMasks, const ElectronicsDelay& electronicsDelay) : mFeeIds(feeIds), mCrateMasks(crateMasks), mElectronicsDelay(electronicsDelay) {} void init(o2::framework::InitContext& ic) { diff --git a/Detectors/TOF/workflow/src/TOFIntegrateClusterSpec.cxx b/Detectors/TOF/workflow/src/TOFIntegrateClusterSpec.cxx index 5e58f5715f63a..2b647f9236611 100644 --- a/Detectors/TOF/workflow/src/TOFIntegrateClusterSpec.cxx +++ b/Detectors/TOF/workflow/src/TOFIntegrateClusterSpec.cxx @@ -174,10 +174,10 @@ o2::framework::DataProcessorSpec getTOFIntegrateClusterSpec(const bool disableWr inputs); std::vector outputs; - outputs.emplace_back(o2::header::gDataOriginTOF, "ITOFCN", 0, Lifetime::Sporadic); - outputs.emplace_back(o2::header::gDataOriginTOF, "ITOFCQ", 0, Lifetime::Sporadic); + outputs.emplace_back(o2::header::gDataOriginTOF, "ITOFCN", 0, Lifetime::Timeframe); + outputs.emplace_back(o2::header::gDataOriginTOF, "ITOFCQ", 0, Lifetime::Timeframe); if (!disableWriter) { - outputs.emplace_back(o2::header::gDataOriginTOF, "ITOFTFID", 0, Lifetime::Sporadic); + outputs.emplace_back(o2::header::gDataOriginTOF, "ITOFTFID", 0, Lifetime::Timeframe); } return DataProcessorSpec{ diff --git a/Detectors/TOF/workflowIO/src/TOFIntegrateClusterWriterSpec.cxx b/Detectors/TOF/workflowIO/src/TOFIntegrateClusterWriterSpec.cxx index 991992783de3a..e230198d668b2 100644 --- a/Detectors/TOF/workflowIO/src/TOFIntegrateClusterWriterSpec.cxx +++ b/Detectors/TOF/workflowIO/src/TOFIntegrateClusterWriterSpec.cxx @@ -30,9 +30,9 @@ DataProcessorSpec getTOFIntegrateClusterWriterSpec() return MakeRootTreeWriterSpec("tof-currents-writer", "o2currents_tof.root", "itofc", - BranchDefinition>{InputSpec{"itofcn", o2::header::gDataOriginTOF, "ITOFCN", 0}, "ITOFCN", 1}, - BranchDefinition>{InputSpec{"itofcq", o2::header::gDataOriginTOF, "ITOFCQ", 0}, "ITOFCQ", 1}, - BranchDefinition{InputSpec{"itoftfid", o2::header::gDataOriginTOF, "ITOFTFID", 0}, "tfID", 1})(); + BranchDefinition>{InputSpec{"itofcn", o2::header::gDataOriginTOF, "ITOFCN", 0, Lifetime::Timeframe}, "ITOFCN", 1}, + BranchDefinition>{InputSpec{"itofcq", o2::header::gDataOriginTOF, "ITOFCQ", 0, Lifetime::Timeframe}, "ITOFCQ", 1}, + BranchDefinition{InputSpec{"itoftfid", o2::header::gDataOriginTOF, "ITOFTFID", 0, Lifetime::Timeframe}, "tfID", 1})(); } } // namespace tof diff --git a/Detectors/TPC/base/include/TPCBase/CDBTypes.h b/Detectors/TPC/base/include/TPCBase/CDBTypes.h index c17d21c8587ae..78ffacbd552fa 100644 --- a/Detectors/TPC/base/include/TPCBase/CDBTypes.h +++ b/Detectors/TPC/base/include/TPCBase/CDBTypes.h @@ -71,6 +71,8 @@ enum class CDBType { /// 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) /// @@ -78,6 +80,9 @@ enum class CDBType { CalScaler, ///< Scaler from IDCs or combined estimator /// 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 }; /// Storage name in CCDB for each calibration and parameter type @@ -129,6 +134,9 @@ const std::unordered_map CDBTypeMap{ // correction maps {CDBType::CalCorrMap, "TPC/Calib/CorrectionMapV2"}, {CDBType::CalCorrMapRef, "TPC/Calib/CorrectionMapRefV2"}, + // correction maps for MC + {CDBType::CalCorrMapMC, "TPC/Calib/CorrectionMapMCV2"}, + {CDBType::CalCorrDerivMapMC, "TPC/Calib/CorrectionMapDerivativeMCV2"}, // derivative map correction {CDBType::CalCorrDerivMap, "TPC/Calib/CorrectionMapDerivativeV2"}, // time series @@ -136,6 +144,9 @@ const std::unordered_map CDBTypeMap{ {CDBType::CalScaler, "TPC/Calib/Scaler"}, // correction maps loader params {CDBType::CorrMapParam, "TPC/Calib/CorrMapParam"}, + // distortion maps + {CDBType::DistortionMapMC, "TPC/Calib/DistortionMapMC"}, + {CDBType::DistortionMapDerivMC, "TPC/Calib/DistortionMapDerivativeMC"}, }; } // namespace o2::tpc diff --git a/Detectors/TPC/calibration/src/CorrectionMapsLoader.cxx b/Detectors/TPC/calibration/src/CorrectionMapsLoader.cxx index 9a2af80934a7e..a9e8ff7971c1a 100644 --- a/Detectors/TPC/calibration/src/CorrectionMapsLoader.cxx +++ b/Detectors/TPC/calibration/src/CorrectionMapsLoader.cxx @@ -67,11 +67,16 @@ void CorrectionMapsLoader::extractCCDBInputs(ProcessingContext& pc) //________________________________________________________ void CorrectionMapsLoader::requestCCDBInputs(std::vector& inputs, std::vector& options, const CorrectionMapsLoaderGloOpts& gloOpts) { - addInput(inputs, {"tpcCorrMap", "TPC", "CorrMap", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrMap), {}, 1)}); // time-dependent if (gloOpts.lumiMode == 0) { + addInput(inputs, {"tpcCorrMap", "TPC", "CorrMap", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrMap), {}, 1)}); // time-dependent addInput(inputs, {"tpcCorrMapRef", "TPC", "CorrMapRef", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrMapRef), {}, 0)}); // load once } else if (gloOpts.lumiMode == 1) { + addInput(inputs, {"tpcCorrMap", "TPC", "CorrMap", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrMap), {}, 1)}); // time-dependent addInput(inputs, {"tpcCorrMapRef", "TPC", "CorrMapRef", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrDerivMap), {}, 1)}); // time-dependent + } else if (gloOpts.lumiMode == 2) { + // for MC corrections + addInput(inputs, {"tpcCorrMap", "TPC", "CorrMap", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrMapMC), {}, 1)}); // time-dependent + addInput(inputs, {"tpcCorrMapRef", "TPC", "CorrMapRef", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrDerivMapMC), {}, 1)}); // time-dependent } else { LOG(fatal) << "Correction mode unknown! Choose either 0 (default) or 1 (derivative map) for flag corrmap-lumi-mode."; } @@ -99,7 +104,7 @@ void CorrectionMapsLoader::addGlobalOptions(std::vector& option { // these are options which should be added at the workflow level, since they modify the inputs of the devices addOption(options, ConfigParamSpec{"lumi-type", o2::framework::VariantType::Int, 0, {"1 = require CTP lumi for TPC correction scaling, 2 = require TPC scalers for TPC correction scaling"}}); - addOption(options, ConfigParamSpec{"corrmap-lumi-mode", o2::framework::VariantType::Int, 0, {"scaling mode: (default) 0 = static + scale * full; 1 = full + scale * derivative"}}); + addOption(options, ConfigParamSpec{"corrmap-lumi-mode", o2::framework::VariantType::Int, 0, {"scaling mode: (default) 0 = static + scale * full; 1 = full + scale * derivative; 2 = full + scale * derivative (for MC)"}}); } //________________________________________________________ diff --git a/Detectors/TPC/calibration/src/VDriftHelper.cxx b/Detectors/TPC/calibration/src/VDriftHelper.cxx index 0df353adba91f..f114084b215b8 100644 --- a/Detectors/TPC/calibration/src/VDriftHelper.cxx +++ b/Detectors/TPC/calibration/src/VDriftHelper.cxx @@ -113,7 +113,7 @@ void VDriftHelper::accountDriftCorrectionITSTPCTgl(const VDriftCorrFact* calib) } if (!mForceParamOffset && mVDTPCITSTgl.refTimeOffset != prevRefTOffs) { // we want to keep the same reference over the run, this should not happen! LOGP(warn, "VDriftHelper: renorming updated TPC refTimeOffset={}/correction={} previous refTimeOffset {}, source: {}", mVDTPCITSTgl.refTimeOffset, mVDTPCITSTgl.timeOffsetCorr, prevRefTOffs, getSourceName()); - mVDTPCITSTgl.normalize(prevRefTOffs); + mVDTPCITSTgl.normalizeOffset(prevRefTOffs); } } } diff --git a/Detectors/TPC/qc/include/TPCQC/Tracks.h b/Detectors/TPC/qc/include/TPCQC/Tracks.h index 509855673ee48..419511753b6cf 100644 --- a/Detectors/TPC/qc/include/TPCQC/Tracks.h +++ b/Detectors/TPC/qc/include/TPCQC/Tracks.h @@ -69,11 +69,13 @@ class Tracks // To set the elementary track cuts void setTrackCuts(float AbsEta = 1., - int nClusterCut = 60, float dEdxTot = 20) + int nClusterCut = 60, float dEdxTot = 20, float cutPtForDCAr = 1.5, float samplingFractionDCAr = 0.1) { mCutAbsEta = AbsEta; mCutMinnCls = nClusterCut; mCutMindEdxTot = dEdxTot; + mCutMinPtDCAr = cutPtForDCAr; + mSamplingFractionDCAr = samplingFractionDCAr; } // Just for backward compatibility with crrent QC, temporary, will be removed in the next PR @@ -96,9 +98,11 @@ class Tracks const std::unordered_map>& getMapHist() const { return mMapHist; } private: - float mCutAbsEta = 1.f; // Eta cut - int mCutMinnCls = 60; // minimum N clusters - float mCutMindEdxTot = 20.f; // dEdxTot min value + float mCutAbsEta = 1.f; // Eta cut + int mCutMinnCls = 60; // minimum N clusters + float mCutMindEdxTot = 20.f; // dEdxTot min value + float mCutMinPtDCAr = 1.5f; // minimum pT for DCAr plots DCAr vs. phi, eta, nCluster + float mSamplingFractionDCAr = 0.1f; // sampling rate for calculation of DCAr std::unordered_map> mMapHist; std::vector mHist1D{}; ///< Initialize vector of 1D histograms std::vector mHist2D{}; ///< Initialize vector of 2D histograms diff --git a/Detectors/TPC/qc/src/Tracks.cxx b/Detectors/TPC/qc/src/Tracks.cxx index 878212c1e0429..86f13d60a00fa 100644 --- a/Detectors/TPC/qc/src/Tracks.cxx +++ b/Detectors/TPC/qc/src/Tracks.cxx @@ -16,6 +16,7 @@ // root includes #include "TFile.h" +#include "TRandom3.h" // o2 includes #include "DataFormatsTPC/TrackTPC.h" @@ -85,6 +86,13 @@ void Tracks::initializeHistograms() for (const auto type : types) { mMapHist[fmt::format("hDCAr_{}", type).data()] = std::make_unique(fmt::format("hDCAr_{}", type).data(), fmt::format("DCAr {};phi;DCAr (cm)", type).data(), 360, 0, o2::math_utils::twoPid(), 250, -10., 10.); } + // DCA vs variables Histograms + mMapHist["hDCArVsPtPos"] = std::make_unique("hDCArVsPtPos", "DCAr Pos;#it{p}_{T}T (GeV/#it{c});DCAr (cm)", logPtBinning.size() - 1, logPtBinning.data(), 250, -10., 10.); + mMapHist["hDCArVsEtaPos"] = std::make_unique("hDCArVsEtaPos", "DCAr Pos;#eta;DCAr (cm)", 400, -2., 2., 250, -10., 10.); + mMapHist["hDCArVsNClsPos"] = std::make_unique("hDCArVsNClsPos", "DCAr Pos;NClusters;DCAr (cm)", 400, -0.5, 399.5, 250, -10., 10.); + mMapHist["hDCArVsPtNeg"] = std::make_unique("hDCArVsPtNeg", "DCAr Neg;#it{p}_{T}T (GeV/#it{c});DCAr (cm)", logPtBinning.size() - 1, logPtBinning.data(), 250, -10., 10.); + mMapHist["hDCArVsEtaNeg"] = std::make_unique("hDCArVsEtaNeg", "DCAr Neg;#eta;DCAr (cm)", 400, -2., 2., 250, -10., 10.); + mMapHist["hDCArVsNClsNeg"] = std::make_unique("hDCArVsNClsNeg", "DCAr Neg;NClusters;DCAr (cm)", 400, -0.5, 399.5, 250, -10., 10.); } //______________________________________________________________________________ void Tracks::resetHistograms() @@ -130,23 +138,46 @@ bool Tracks::processTrack(const o2::tpc::TrackTPC& track) auto propagator = o2::base::Propagator::Instance(true); const int type = (track.getQ2Pt() < 0) + 2 * track.hasCSideClustersOnly(); auto dcaHist = mMapHist[fmt::format("hDCAr_{}", types[type]).data()].get(); + const std::string signType((sign < 0) ? "Neg" : "Pos"); + auto dcaHistPT = mMapHist["hDCArVsPt" + signType].get(); + auto dcaHistEta = mMapHist["hDCArVsEta" + signType].get(); + auto dcaHistNCluster = mMapHist["hDCArVsNCls" + signType].get(); - if (propagator->getMatLUT() && propagator->hasMagFieldSet()) { - // ---| fill DCA histos |--- - o2::gpu::gpustd::array dca; - const o2::math_utils::Point3D refPoint{0, 0, 0}; - o2::track::TrackPar propTrack(track); - if (propagator->propagateToDCABxByBz(refPoint, propTrack, 2.f, o2::base::Propagator::MatCorrType::USEMatCorrLUT, &dca)) { - const auto phi = o2::math_utils::to02PiGen(track.getPhi()); - dcaHist->Fill(phi, dca[0]); - } - } else { - static bool reported = false; - if (!reported) { - LOGP(error, "o2::base::Propagator not properly initialized, MatLUT ({}) and / or Field ({}) missing, will not fill DCA histograms", (void*)propagator->getMatLUT(), (void*)propagator->hasMagFieldSet()); - reported = true; + // set-up sampling for the DCA calculation + Double_t sampleProb = 2; + + if (mSamplingFractionDCAr > 0) { // for now no SEED is given. + TRandom3 randomGenerator(0); + sampleProb = randomGenerator.Uniform(1); + } + + if (sampleProb > (Double_t)(1. - mSamplingFractionDCAr)) { + + if (propagator->getMatLUT() && propagator->hasMagFieldSet()) { + // ---| fill DCA histos |--- + o2::gpu::gpustd::array dca; + const o2::math_utils::Point3D refPoint{0, 0, 0}; + o2::track::TrackPar propTrack(track); + if (propagator->propagateToDCABxByBz(refPoint, propTrack, 2.f, o2::base::Propagator::MatCorrType::USEMatCorrLUT, &dca)) { + const auto phi = o2::math_utils::to02PiGen(track.getPhi()); + dcaHistPT->Fill(pt, dca[0]); + if (pt > mCutMinPtDCAr) { + dcaHist->Fill(phi, dca[0]); + dcaHistEta->Fill(eta, dca[0]); + dcaHistNCluster->Fill(nCls, dca[0]); + } + } else { + static bool reported = false; + if (!reported) { + LOGP(error, "o2::base::Propagator not properly initialized, MatLUT ({}) and / or Field ({}) missing, will not fill DCA histograms", (void*)propagator->getMatLUT(), (void*)propagator->hasMagFieldSet()); + reported = true; + } + dcaHist->SetTitle(fmt::format("DCAr {} o2::base::Propagator not properly initialized", types[type]).data()); + dcaHistPT->SetTitle(fmt::format("DCAr p_{T} {} o2::base::Propagator not properly initialized", signType).data()); + dcaHistEta->SetTitle(fmt::format("DCAr eta {} o2::base::Propagator not properly initialized", signType).data()); + dcaHistNCluster->SetTitle(fmt::format("DCAr nClusters {} o2::base::Propagator not properly initialized", signType).data()); + } } - dcaHist->SetTitle(fmt::format("DCAr {} o2::base::Propagator not properly initialized", types[type]).data()); } if (hasASideOnly == 1) { diff --git a/Detectors/TPC/simulation/CMakeLists.txt b/Detectors/TPC/simulation/CMakeLists.txt index 1f433469ab769..510d42d2d6d85 100644 --- a/Detectors/TPC/simulation/CMakeLists.txt +++ b/Detectors/TPC/simulation/CMakeLists.txt @@ -22,7 +22,7 @@ o2_add_library(TPCSimulation src/SAMPAProcessing.cxx src/IDCSim.cxx PUBLIC_LINK_LIBRARIES O2::DetectorsBase O2::SimulationDataFormat - O2::TPCBase O2::TPCSpaceCharge + O2::TPCBase O2::TPCSpaceCharge O2::TPCCalibration ROOT::Physics) o2_target_root_dictionary(TPCSimulation diff --git a/Detectors/TPC/simulation/include/TPCSimulation/Digitizer.h b/Detectors/TPC/simulation/include/TPCSimulation/Digitizer.h index 359dc5d64b72d..e2c7e9ec7d175 100644 --- a/Detectors/TPC/simulation/include/TPCSimulation/Digitizer.h +++ b/Detectors/TPC/simulation/include/TPCSimulation/Digitizer.h @@ -53,7 +53,7 @@ enum class SCDistortionType : int; class Digitizer { public: - using SC = SpaceCharge; + using SC = SpaceCharge; /// Default constructor Digitizer(); @@ -121,6 +121,9 @@ class Digitizer /// \param spaceCharge unique pointer to spaceCharge object void setUseSCDistortions(SC* spaceCharge); + /// \param spaceCharge unique pointer to spaceCharge object + void setSCDistortionsDerivative(SC* spaceCharge); + /// Enable the use of space-charge distortions by providing global distortions and global corrections stored in a ROOT file /// The storage of the values should be done by the methods provided in the SpaceCharge class /// \param file containing distortions @@ -129,17 +132,26 @@ class Digitizer void setVDrift(float v) { mVDrift = v; } void setTDriftOffset(float t) { mTDriftOffset = t; } + void setDistortionScaleType(int distortionScaleType) { mDistortionScaleType = distortionScaleType; } + int getDistortionScaleType() const { return mDistortionScaleType; } + void setLumiScaleFactor(); + void setMeanLumiDistortions(float meanLumi); + void setMeanLumiDistortionsDerivative(float meanLumi); + private: - DigitContainer mDigitContainer; ///< Container for the Digits - std::unique_ptr mSpaceCharge; ///< Handler of space-charge distortions - Sector mSector = -1; ///< ID of the currently processed sector - double mEventTime = 0.f; ///< Time of the currently processed event - double mOutputDigitTimeOffset = 0; ///< Time of the first IR sampled in the digitizer - float mVDrift = 0; ///< VDrift for current timestamp - float mTDriftOffset = 0; ///< drift time additive offset in \mus - bool mIsContinuous; ///< Switch for continuous readout - bool mUseSCDistortions = false; ///< Flag to switch on the use of space-charge distortions - ClassDefNV(Digitizer, 1); + DigitContainer mDigitContainer; ///< Container for the Digits + std::unique_ptr mSpaceCharge; ///< Handler of full distortions (static + IR dependant) + std::unique_ptr mSpaceChargeDer; ///< Handler of reference static distortions + Sector mSector = -1; ///< ID of the currently processed sector + double mEventTime = 0.f; ///< Time of the currently processed event + double mOutputDigitTimeOffset = 0; ///< Time of the first IR sampled in the digitizer + float mVDrift = 0; ///< VDrift for current timestamp + float mTDriftOffset = 0; ///< drift time additive offset in \mus + bool mIsContinuous; ///< Switch for continuous readout + bool mUseSCDistortions = false; ///< Flag to switch on the use of space-charge distortions + int mDistortionScaleType = 0; ///< type=0: no scaling of distortions, type=1 distortions without any scaling, type=2 distortions scaling with lumi + float mLumiScaleFactor = 0; ///< value used to scale the derivative map + ClassDefNV(Digitizer, 2); }; } // namespace tpc } // namespace o2 diff --git a/Detectors/TPC/simulation/src/Digitizer.cxx b/Detectors/TPC/simulation/src/Digitizer.cxx index 7bac9e6ebb3e0..38968b3b38bb9 100644 --- a/Detectors/TPC/simulation/src/Digitizer.cxx +++ b/Detectors/TPC/simulation/src/Digitizer.cxx @@ -27,6 +27,7 @@ #include "TPCBase/CDBInterface.h" #include "TPCSpaceCharge/SpaceCharge.h" #include "TPCBase/Mapper.h" +#include "TPCCalibration/CorrMapParam.h" #include @@ -40,10 +41,6 @@ Digitizer::Digitizer() = default; void Digitizer::init() { - // Calculate distortion lookup tables if initial space-charge density is provided - if (mUseSCDistortions) { - mSpaceCharge->init(); - } auto& gemAmplification = GEMAmplification::instance(); gemAmplification.updateParameters(); auto& electronTransport = ElectronTransport::instance(); @@ -83,8 +80,10 @@ void Digitizer::process(const std::vector& hits, GlobalPosition3D posEle(eh.GetX(), eh.GetY(), eh.GetZ()); // Distort the electron position in case space-charge distortions are used - if (mUseSCDistortions) { + if (mDistortionScaleType == 1) { mSpaceCharge->distortElectron(posEle); + } else if (mDistortionScaleType == 2) { + mSpaceCharge->distortElectron(posEle, mSpaceChargeDer.get(), mLumiScaleFactor); } /// Remove electrons that end up more than three sigma of the hit's average diffusion away from the current sector @@ -190,6 +189,15 @@ void Digitizer::setUseSCDistortions(SC* spaceCharge) { mUseSCDistortions = true; mSpaceCharge.reset(spaceCharge); + mSpaceCharge->initAfterReadingFromFile(); + mSpaceCharge->printMetaData(); +} + +void Digitizer::setSCDistortionsDerivative(SC* spaceCharge) +{ + mSpaceChargeDer.reset(spaceCharge); + mSpaceChargeDer->initAfterReadingFromFile(); + mSpaceChargeDer->printMetaData(); } void Digitizer::setUseSCDistortions(std::string_view finp) @@ -213,3 +221,19 @@ void Digitizer::setStartTime(double time) sampaProcessing.updateParameters(mVDrift); mDigitContainer.setStartTime(sampaProcessing.getTimeBinFromTime(time - mOutputDigitTimeOffset)); } + +void Digitizer::setLumiScaleFactor() +{ + mLumiScaleFactor = (CorrMapParam::Instance().lumiInst - mSpaceCharge->getMeanLumi()) / mSpaceChargeDer->getMeanLumi(); + LOGP(info, "Setting Lumi scale factor: lumiInst: {} lumi mean: {} lumi mean derivative: {} lumi scale factor: {}", CorrMapParam::Instance().lumiInst, mSpaceCharge->getMeanLumi(), mSpaceChargeDer->getMeanLumi(), mLumiScaleFactor); +} + +void Digitizer::setMeanLumiDistortions(float meanLumi) +{ + mSpaceCharge->setMeanLumi(meanLumi); +} + +void Digitizer::setMeanLumiDistortionsDerivative(float meanLumi) +{ + mSpaceChargeDer->setMeanLumi(meanLumi); +} diff --git a/Detectors/TPC/spacecharge/include/TPCSpaceCharge/SpaceCharge.h b/Detectors/TPC/spacecharge/include/TPCSpaceCharge/SpaceCharge.h index 014be1f2cc50e..dc6bb1be3aba9 100644 --- a/Detectors/TPC/spacecharge/include/TPCSpaceCharge/SpaceCharge.h +++ b/Detectors/TPC/spacecharge/include/TPCSpaceCharge/SpaceCharge.h @@ -572,7 +572,9 @@ class SpaceCharge /// Distort electron position using distortion lookup tables /// \param point 3D coordinates of the electron - void distortElectron(GlobalPosition3D& point) const; + /// \param scSCale other sc object which is used for scaling of the distortions + /// \param scale scaling value + void distortElectron(GlobalPosition3D& point, const SpaceCharge* scSCale = nullptr, float scale = 0) const; /// set the distortions directly from a look up table /// \param distdZ distortions in z direction @@ -1170,12 +1172,24 @@ class SpaceCharge /// substract global corrections from other sc object (global corrections -= other.global corrections) /// can be used to calculate the derivative: (this - other)/normalization /// for normalization see scaleCorrections() - void substractGlobalCorrections(const SpaceCharge& otherSC, const Side side); + void subtractGlobalCorrections(const SpaceCharge& otherSC, const Side side); + + /// substract global distortions from other sc object (global distortions -= other.global distortions) + /// can be used to calculate the derivative: (this - other)/normalization + void subtractGlobalDistortions(const SpaceCharge& otherSC, const Side side); /// scale corrections by factor /// \param scaleFac global corrections are multiplied by this factor void scaleCorrections(const float scaleFac, const Side side); + /// setting meta data for this object + void setMetaData(const SCMetaData& meta) { mMeta = meta; } + const auto& getMetaData() const { return mMeta; } + void printMetaData() const { mMeta.print(); } + float getMeanLumi() const { return mMeta.meanLumi; } + void setMeanLumi(float lumi) { mMeta.meanLumi = lumi; } + void initAfterReadingFromFile(); + private: ParamSpaceCharge mParamGrid{}; ///< parameters of the grid on which the calculations are performed inline static int sNThreads{getOMPMaxThreads()}; /// mAnaDistCorr; ///< analytical distortions and corrections bool mUseAnaDistCorr{false}; ///< flag if analytical distortions will be used in the distortElectron() and getCorrections() function BField mBField{}; /// collisionTypes{"PP", "Pb-Pb"}; + const std::array slumiSource{"CTP", "IDC"}; + if (collisionType < collisionTypes.size()) { + LOGP(info, "meanLumi: {}, source of lumi: {}, IR: {}kHz, run: {}, collisionType {}", meanLumi, slumiSource[lumiSource], ir, run, collisionTypes[collisionType]); + } else { + LOGP(info, "Specified collision type {} not allowed", collisionType); + } + } + + float meanLumi = 0; ///< mean lumi the sc object corresponds to + float ir = 0; ///< IR + int run = 0; ///< run number this object anchored to to + int collisionType = 0; ///< 0=PP, 1-Pb-Pb + LumiType lumiSource{}; ///< source of luminosity + + private: + ClassDefNV(SCMetaData, 1); +}; + /// /// this class contains an analytical description of the space charge, potential and the electric fields. /// The analytical functions can be used to test the poisson solver and the caluclation of distortions/corrections. diff --git a/Detectors/TPC/spacecharge/src/SpaceCharge.cxx b/Detectors/TPC/spacecharge/src/SpaceCharge.cxx index 3ed00345159bc..d3fa62b94a1c9 100644 --- a/Detectors/TPC/spacecharge/src/SpaceCharge.cxx +++ b/Detectors/TPC/spacecharge/src/SpaceCharge.cxx @@ -19,7 +19,6 @@ #include "fmt/core.h" #include "Framework/Logger.h" #include "TPCSpaceCharge/PoissonSolver.h" -#include "Framework/Logger.h" #include "TGeoGlobalMagField.h" #include "TPCBase/ParameterGas.h" #include "TPCBase/ParameterElectronics.h" @@ -709,7 +708,7 @@ void SpaceCharge::calcGlobalDistWithGlobalCorrIterative(const DistCorrInt const DataT phi = getPhiVertex(iPhi, side); for (unsigned int iR = 0; iR < mParamGrid.NRVertices; ++iR) { const DataT radius = getRVertex(iR, side); - for (unsigned int iZ = 0; iZ < mParamGrid.NZVertices; ++iZ) { + for (unsigned int iZ = 1; iZ < mParamGrid.NZVertices; ++iZ) { const DataT z = getZVertex(iZ, side); unsigned int nearestiZ = iZ; @@ -789,6 +788,11 @@ void SpaceCharge::calcGlobalDistWithGlobalCorrIterative(const DistCorrInt mGlobalDistdZ[side](iZ, iR, iPhi) = -corrdZ; } } + for (unsigned int iR = 0; iR < mParamGrid.NRVertices; ++iR) { + mGlobalDistdR[side](0, iR, iPhi) = 3 * (mGlobalDistdR[side](1, iR, iPhi) - mGlobalDistdR[side](2, iR, iPhi)) + mGlobalDistdR[side](3, iR, iPhi); + mGlobalDistdRPhi[side](0, iR, iPhi) = 3 * (mGlobalDistdRPhi[side](1, iR, iPhi) - mGlobalDistdRPhi[side](2, iR, iPhi)) + mGlobalDistdRPhi[side](3, iR, iPhi); + mGlobalDistdZ[side](0, iR, iPhi) = 3 * (mGlobalDistdZ[side](1, iR, iPhi) - mGlobalDistdZ[side](2, iR, iPhi)) + mGlobalDistdZ[side](3, iR, iPhi); + } } } @@ -1516,7 +1520,7 @@ void SpaceCharge::correctElectron(GlobalPosition3D& point) } template -void SpaceCharge::distortElectron(GlobalPosition3D& point) const +void SpaceCharge::distortElectron(GlobalPosition3D& point, const SpaceCharge* scSCale, float scale) const { DataT distX{}; DataT distY{}; @@ -1525,6 +1529,18 @@ void SpaceCharge::distortElectron(GlobalPosition3D& point) const // get the distortions for input coordinate getDistortions(point.X(), point.Y(), point.Z(), side, distX, distY, distZ); + DataT distXTmp{}; + DataT distYTmp{}; + DataT distZTmp{}; + + // scale distortions if requested + if (scSCale && scale != 0) { + scSCale->getDistortions(point.X(), point.Y(), point.Z(), side, distXTmp, distYTmp, distZTmp); + distX += distXTmp * scale; + distY += distYTmp * scale; + distZ += distZTmp * scale; + } + GPUCA_DEBUG_STREAMER_CHECK(if (o2::utils::DebugStreamer::checkStream(o2::utils::StreamFlags::streamDistortionsSC)) { GlobalPosition3D pos(point); float phi = std::atan2(pos.Y(), pos.X()); @@ -1543,6 +1559,10 @@ void SpaceCharge::distortElectron(GlobalPosition3D& point) const << "distX=" << distX << "distY=" << distY << "distZ=" << distZ + << "distXDer=" << distXTmp + << "distYDer=" << distYTmp + << "distZDer=" << distZTmp + << "scale=" << scale << "\n"; }) @@ -1734,8 +1754,10 @@ void SpaceCharge::getDistortionsCyl(const std::vector& z, const st template void SpaceCharge::getDistortions(const DataT x, const DataT y, const DataT z, const Side side, DataT& distX, DataT& distY, DataT& distZ) const { + DataT zClamped = regulateZ(z, side); + if (mUseAnaDistCorr) { - getDistortionsAnalytical(x, y, z, side, distX, distY, distZ); + getDistortionsAnalytical(x, y, zClamped, side, distX, distY, distZ); } else { // convert cartesian to polar const DataT radius = getRadiusFromCartesian(x, y); @@ -1743,11 +1765,12 @@ void SpaceCharge::getDistortions(const DataT x, const DataT y, const Data DataT distR{}; DataT distRPhi{}; - getDistortionsCyl(z, radius, phi, side, distZ, distR, distRPhi); + DataT rClamped = regulateR(radius, side); + getDistortionsCyl(zClamped, rClamped, phi, side, distZ, distR, distRPhi); // Calculate distorted position - const DataT radiusDist = radius + distR; - const DataT phiDist = phi + distRPhi / radius; + const DataT radiusDist = rClamped + distR; + const DataT phiDist = phi + distRPhi / rClamped; distX = getXFromPolar(radiusDist, phiDist) - x; // difference between distorted and original x coordinate distY = getYFromPolar(radiusDist, phiDist) - y; // difference between distorted and original y coordinate @@ -2268,7 +2291,13 @@ void SpaceCharge::dumpToTree(const char* outFileName, const Side side, co DataT corrZ{}; DataT corrR{}; DataT corrRPhi{}; - getCorrectionsCyl(zPos, rPos, phiPos, side, corrZ, corrR, corrRPhi); + // getCorrectionsCyl(zPos, rPos, phiPos, side, corrZ, corrR, corrRPhi); + + const DataT zDistorted = zPos + distZ; + const DataT radiusDistorted = rPos + distR; + const DataT phiDistorted = regulatePhi(phiPos + distRPhi / rPos, side); + getCorrectionsCyl(zDistorted, radiusDistorted, phiDistorted, side, corrZ, corrR, corrRPhi); + corrRPhi *= rPos / radiusDistorted; DataT lcorrZ{}; DataT lcorrR{}; @@ -3046,12 +3075,13 @@ void SpaceCharge::dumpMetaData(std::string_view file, std::string_view op dfStore = dfStore.DefineSlotEntry("grid_A", [&helperA = helperA](unsigned int, ULong64_t entry) { return helperA; }); dfStore = dfStore.DefineSlotEntry("grid_C", [&helperC = helperC](unsigned int, ULong64_t entry) { return helperC; }); dfStore = dfStore.DefineSlotEntry("BField", [field = mBField.getBField()](unsigned int, ULong64_t entry) { return field; }); + dfStore = dfStore.DefineSlotEntry("metaInf", [meta = mMeta](unsigned int, ULong64_t entry) { return meta; }); // write to TTree ROOT::RDF::RSnapshotOptions opt; opt.fMode = option; opt.fOverwriteIfExists = true; // overwrite if already exists - dfStore.Snapshot("meta", file, {"paramsC", "grid_A", "grid_C", "BField"}, opt); + dfStore.Snapshot("meta", file, {"paramsC", "grid_A", "grid_C", "BField", "metaInf"}, opt); } template @@ -3079,6 +3109,15 @@ void SpaceCharge::readMetaData(std::string_view file) ROOT::RDataFrame dFrame("meta", file); dFrame.Foreach(readMeta, {"paramsC", "grid_A", "grid_C", "BField"}); + + const auto& cols = dFrame.GetColumnNames(); + if (std::find(cols.begin(), cols.end(), "metaInf") != cols.end()) { + auto readMetaInf = [&mMeta = mMeta](const SCMetaData& meta) { + mMeta = meta; + }; + dFrame.Foreach(readMetaInf, {"metaInf"}); + } + LOGP(info, "Setting meta data: mC0={} mC1={} mC2={}", mC0, mC1, mC2); mReadMetaData = true; } @@ -3608,13 +3647,21 @@ void SpaceCharge::fillROCMisalignment(const std::vector& indicesT } template -void SpaceCharge::substractGlobalCorrections(const SpaceCharge& otherSC, const Side side) +void SpaceCharge::subtractGlobalCorrections(const SpaceCharge& otherSC, const Side side) { mGlobalCorrdR[side] -= otherSC.mGlobalCorrdR[side]; mGlobalCorrdZ[side] -= otherSC.mGlobalCorrdZ[side]; mGlobalCorrdRPhi[side] -= otherSC.mGlobalCorrdRPhi[side]; } +template +void SpaceCharge::subtractGlobalDistortions(const SpaceCharge& otherSC, const Side side) +{ + mGlobalDistdR[side] -= otherSC.mGlobalDistdR[side]; + mGlobalDistdZ[side] -= otherSC.mGlobalDistdZ[side]; + mGlobalDistdRPhi[side] -= otherSC.mGlobalDistdRPhi[side]; +} + template void SpaceCharge::scaleCorrections(const float val, const Side side) { @@ -3694,6 +3741,13 @@ void SpaceCharge::scaleChargeDensityStack(const float scalingFactor, cons } } +template +void SpaceCharge::initAfterReadingFromFile() +{ + mGrid3D[Side::A] = RegularGrid(GridProp::ZMIN, GridProp::RMIN, GridProp::PHIMIN, getSign(Side::A) * GridProp::getGridSpacingZ(mParamGrid.NZVertices), GridProp::getGridSpacingR(mParamGrid.NRVertices), GridProp::getGridSpacingPhi(mParamGrid.NPhiVertices), mParamGrid); + mGrid3D[Side::C] = RegularGrid(GridProp::ZMIN, GridProp::RMIN, GridProp::PHIMIN, getSign(Side::C) * GridProp::getGridSpacingZ(mParamGrid.NZVertices), GridProp::getGridSpacingR(mParamGrid.NRVertices), GridProp::getGridSpacingPhi(mParamGrid.NPhiVertices), mParamGrid); +} + using DataTD = double; template class o2::tpc::SpaceCharge; diff --git a/Detectors/TPC/spacecharge/src/TPCSpacechargeLinkDef.h b/Detectors/TPC/spacecharge/src/TPCSpacechargeLinkDef.h index 7574b6896c315..b1bb33dc1aeb7 100644 --- a/Detectors/TPC/spacecharge/src/TPCSpacechargeLinkDef.h +++ b/Detectors/TPC/spacecharge/src/TPCSpacechargeLinkDef.h @@ -47,4 +47,5 @@ #pragma link C++ class o2::tpc::AnalyticalDistCorr < double> + ; #pragma link C++ class o2::tpc::AnalyticalDistCorr < float> + ; #pragma link C++ struct o2::tpc::ParamSpaceCharge + ; +#pragma link C++ struct o2::tpc::SCMetaData + ; #endif diff --git a/Detectors/TPC/workflow/src/TPCIntegrateClusterSpec.cxx b/Detectors/TPC/workflow/src/TPCIntegrateClusterSpec.cxx index f06782af16b69..46be10dbc2c69 100644 --- a/Detectors/TPC/workflow/src/TPCIntegrateClusterSpec.cxx +++ b/Detectors/TPC/workflow/src/TPCIntegrateClusterSpec.cxx @@ -185,9 +185,9 @@ o2::framework::DataProcessorSpec getTPCIntegrateClusterSpec(const bool disableWr inputs); std::vector outputs; - outputs.emplace_back(o2::header::gDataOriginTPC, getDataDescriptionTPCC(), 0, Lifetime::Sporadic); + outputs.emplace_back(o2::header::gDataOriginTPC, getDataDescriptionTPCC(), 0, Lifetime::Timeframe); if (!disableWriter) { - outputs.emplace_back(o2::header::gDataOriginTPC, getDataDescriptionTPCTFId(), 0, Lifetime::Sporadic); + outputs.emplace_back(o2::header::gDataOriginTPC, getDataDescriptionTPCTFId(), 0, Lifetime::Timeframe); } return DataProcessorSpec{ diff --git a/Detectors/TPC/workflow/src/TPCIntegrateClusterWriterSpec.cxx b/Detectors/TPC/workflow/src/TPCIntegrateClusterWriterSpec.cxx index 838efc6b8dc3f..40a5ec3e5cc06 100644 --- a/Detectors/TPC/workflow/src/TPCIntegrateClusterWriterSpec.cxx +++ b/Detectors/TPC/workflow/src/TPCIntegrateClusterWriterSpec.cxx @@ -33,8 +33,8 @@ DataProcessorSpec getTPCIntegrateClusterWriterSpec() return MakeRootTreeWriterSpec("tpc-currents-writer", "o2currents_tpc.root", "itpcc", - BranchDefinition{InputSpec{"itpcc", o2::header::gDataOriginTPC, getDataDescriptionTPCC(), 0}, "ITPCC", 1}, - BranchDefinition{InputSpec{"itpctfid", o2::header::gDataOriginTPC, "ITPCTFID", 0}, "tfID", 1})(); + BranchDefinition{InputSpec{"itpcc", o2::header::gDataOriginTPC, getDataDescriptionTPCC(), 0, Lifetime::Sporadic}, "ITPCC", 1}, + BranchDefinition{InputSpec{"itpctfid", o2::header::gDataOriginTPC, "ITPCTFID", 0, Lifetime::Sporadic}, "tfID", 1})(); } } // namespace tpc diff --git a/Detectors/TPC/workflow/src/TPCTimeSeriesSpec.cxx b/Detectors/TPC/workflow/src/TPCTimeSeriesSpec.cxx index 6adb485a45477..4dfa67a7b826b 100644 --- a/Detectors/TPC/workflow/src/TPCTimeSeriesSpec.cxx +++ b/Detectors/TPC/workflow/src/TPCTimeSeriesSpec.cxx @@ -43,6 +43,8 @@ #include "DataFormatsTPC/PIDResponse.h" #include "DataFormatsITS/TrackITS.h" #include "TROOT.h" +#include "ReconstructionDataFormats/MatchInfoTOF.h" +#include "DataFormatsTOF/Cluster.h" using namespace o2::globaltracking; using GTrackID = o2::dataformats::GlobalTrackID; @@ -170,7 +172,10 @@ class TPCTimeSeries : public Task auto primMatchedTracks = mTPCOnly ? gsl::span() : recoData.getPrimaryVertexMatchedTracks(); // Global ID's for associated tracks auto primMatchedTracksRef = mTPCOnly ? gsl::span() : recoData.getPrimaryVertexMatchedTrackRefs(); // references from vertex to these track IDs - LOGP(info, "Processing {} vertices, {} primary matched vertices, {} TPC tracks, {} ITS tracks, {} ITS-TPC tracks", vertices.size(), primMatchedTracks.size(), tracksTPC.size(), tracksITS.size(), tracksITSTPC.size()); + // TOF clusters + const auto& tofClusters = mTPCOnly ? gsl::span() : recoData.getTOFClusters(); + + LOGP(info, "Processing {} vertices, {} primary matched vertices, {} TPC tracks, {} ITS tracks, {} ITS-TPC tracks, {} TOF clusters", vertices.size(), primMatchedTracks.size(), tracksTPC.size(), tracksITS.size(), tracksITSTPC.size(), tofClusters.size()); // calculate mean vertex, RMS and count vertices auto indicesITSTPC_vtx = processVertices(vertices, primMatchedTracks, primMatchedTracksRef, recoData); @@ -186,6 +191,24 @@ class TPCTimeSeries : public Task indicesITSTPC[tracksITSTPC[i].getRefTPC().getIndex()] = {i, idxVtx}; } + std::vector> idxTPCTrackToTOFCluster; // store for each tpc track index the index to the TOF cluster + + // get matches to TOF in case skimmed data is produced + if (mUnbinnedWriter) { + idxTPCTrackToTOFCluster = std::vector>(tracksTPC.size(), {-1, -999, -999}); + const std::vector> tofMatches{recoData.getTPCTOFMatches(), recoData.getTPCTRDTOFMatches(), recoData.getITSTPCTOFMatches(), recoData.getITSTPCTRDTOFMatches()}; + + // loop over ITS-TPC-TRD-TOF and ITS-TPC-TOF tracks an store for each ITS-TPC track the TOF track index + for (const auto& tofMatch : tofMatches) { + for (const auto& tpctofmatch : tofMatch) { + auto refTPC = recoData.getTPCContributorGID(tpctofmatch.getTrackRef()); + if (refTPC.isIndexSet()) { + idxTPCTrackToTOFCluster[refTPC] = {tpctofmatch.getIdxTOFCl(), tpctofmatch.getDXatTOF(), tpctofmatch.getDZatTOF()}; + } + } + } + } + // find nearest vertex of tracks which have no vertex assigned findNearesVertex(tracksTPC, vertices); @@ -358,7 +381,7 @@ class TPCTimeSeries : public Task auto myThread = [&](int iThread) { for (size_t i = iThread; i < loopEnd; i += mNThreads) { if (acceptTrack(tracksTPC[i])) { - fillDCA(tracksTPC, tracksITSTPC, vertices, i, iThread, indicesITSTPC, tracksITS); + fillDCA(tracksTPC, tracksITSTPC, vertices, i, iThread, indicesITSTPC, tracksITS, idxTPCTrackToTOFCluster, tofClusters); } } }; @@ -375,7 +398,7 @@ class TPCTimeSeries : public Task auto myThread = [&](int iThread) { for (size_t i = iThread; i < loopEnd; i += mNThreads) { if (acceptTrack(tracksTPC[i])) { - fillDCA(tracksTPC, tracksITSTPC, vertices, i, iThread, indicesITSTPC, tracksITS); + fillDCA(tracksTPC, tracksITSTPC, vertices, i, iThread, indicesITSTPC, tracksITS, idxTPCTrackToTOFCluster, tofClusters); } } }; @@ -1001,9 +1024,10 @@ class TPCTimeSeries : public Task return true; } - void fillDCA(const gsl::span tracksTPC, const gsl::span tracksITSTPC, const gsl::span vertices, const int iTrk, const int iThread, const std::unordered_map>& indicesITSTPC, const gsl::span tracksITS) + void fillDCA(const gsl::span tracksTPC, const gsl::span tracksITSTPC, const gsl::span vertices, const int iTrk, const int iThread, const std::unordered_map>& indicesITSTPC, const gsl::span tracksITS, const std::vector>& idxTPCTrackToTOFCluster, const gsl::span tofClusters) { - TrackTPC track = tracksTPC[iTrk]; + o2::track::TrackParCov track = tracksTPC[iTrk]; + const auto& trackFull = tracksTPC[iTrk]; // propagate track to the DCA and fill in slice auto propagator = o2::base::Propagator::Instance(); @@ -1055,11 +1079,11 @@ class TPCTimeSeries : public Task const auto vertex = (idxITSTPC.back() != -1) ? vertices[idxITSTPC.back()] : ((mNearestVtxTPC[iTrk] != -1) ? vertices[mNearestVtxTPC[iTrk]] : o2::dataformats::PrimaryVertex{}); // calculate DCAz: (time TPC track - time vertex) * vDrift + sign_side * vertexZ - const float signSide = track.hasCSideClustersOnly() ? -1 : 1; // invert sign for C-side - const float dcaZFromDeltaTime = (vertex.getTimeStamp().getTimeStamp() == 0) ? 0 : (o2::tpc::ParameterElectronics::Instance().ZbinWidth * track.getTime0() - vertex.getTimeStamp().getTimeStamp()) * mVDrift + signSide * vertex.getZ(); + const float signSide = trackFull.hasCSideClustersOnly() ? -1 : 1; // invert sign for C-side + const float dcaZFromDeltaTime = (vertex.getTimeStamp().getTimeStamp() == 0) ? 0 : (o2::tpc::ParameterElectronics::Instance().ZbinWidth * trackFull.getTime0() - vertex.getTimeStamp().getTimeStamp()) * mVDrift + signSide * vertex.getZ(); // for weight of DCA - const float resCl = std::min(track.getNClusters(), static_cast(Mapper::PADROWS)) / static_cast(Mapper::PADROWS); + const float resCl = std::min(trackFull.getNClusters(), static_cast(Mapper::PADROWS)) / static_cast(Mapper::PADROWS); const float div = (resCl * track.getPt()); if (div == 0) { @@ -1087,15 +1111,15 @@ class TPCTimeSeries : public Task } const float chi2Match = (chi2 > 0) ? std::sqrt(chi2) : -1; - const float sqrtChi2TPC = (track.getChi2() > 0) ? std::sqrt(track.getChi2()) : 0; - const float nClTPC = track.getNClusters(); + const float sqrtChi2TPC = (trackFull.getChi2() > 0) ? std::sqrt(trackFull.getChi2()) : 0; + const float nClTPC = trackFull.getNClusters(); // const float dedx = mUseQMax ? track.getdEdx().dEdxMaxTPC : track.getdEdx().dEdxTotTPC; - const float dedxRatioqTot = (track.getdEdx().dEdxTotTPC > 0) ? (mMIPdEdx / track.getdEdx().dEdxTotTPC) : -1; - const float dedxRatioqMax = (track.getdEdx().dEdxMaxTPC > 0) ? (mMIPdEdx / track.getdEdx().dEdxMaxTPC) : -1; + const float dedxRatioqTot = (trackFull.getdEdx().dEdxTotTPC > 0) ? (mMIPdEdx / trackFull.getdEdx().dEdxTotTPC) : -1; + const float dedxRatioqMax = (trackFull.getdEdx().dEdxMaxTPC > 0) ? (mMIPdEdx / trackFull.getdEdx().dEdxMaxTPC) : -1; - const auto dedxQTotVars = getdEdxVars(0, track); - const auto dedxQMaxVars = getdEdxVars(1, track); + const auto dedxQTotVars = getdEdxVars(0, trackFull); + const auto dedxQMaxVars = getdEdxVars(1, trackFull); // make check to avoid crash in case no or less ITS tracks have been found! const int idxITSTrack = (hasITSTPC && (gID == o2::dataformats::GlobalTrackID::Source::ITS)) ? tracksITSTPC[idxITSTPC.front()].getRefITS().getIndex() : -1; @@ -1108,9 +1132,9 @@ class TPCTimeSeries : public Task } sigmaY2 = track.getSigmaY2(); sigmaZ2 = track.getSigmaZ2(); - if (track.hasCSideClustersOnly()) { + if (trackFull.hasCSideClustersOnly()) { mBufferVals[iThread].front().emplace_back(Side::C, tglBin, phiBin, qPtBin, multBin, dca[0], dcaZFromDeltaTime, dcarW, dedxRatioqTot, dedxRatioqMax, sqrtChi2TPC, nClTPC, gID, chi2Match, hasITSTPC, nClITS, chi2ITS, dedxQTotVars, dedxQMaxVars, sigmaY2, sigmaZ2); - } else if (track.hasASideClustersOnly()) { + } else if (trackFull.hasASideClustersOnly()) { mBufferVals[iThread].front().emplace_back(Side::A, tglBin, phiBin, qPtBin, multBin, dca[0], dcaZFromDeltaTime, dcarW, dedxRatioqTot, dedxRatioqMax, sqrtChi2TPC, nClTPC, gID, chi2Match, hasITSTPC, nClITS, chi2ITS, dedxQTotVars, dedxQMaxVars, sigmaY2, sigmaZ2); } @@ -1150,9 +1174,9 @@ class TPCTimeSeries : public Task dcaITSTPCTmp[1] = -1; } - if (track.hasCSideClustersOnly()) { + if (trackFull.hasCSideClustersOnly()) { mBufferVals[iThread].back().emplace_back_ITSTPC(Side::C, tglBin, phiBin, qPtBin, multBin, dca[0], dcaZFromDeltaTime, dcarW, dedxRatioqTot, dedxRatioqMax, sqrtChi2TPC, nClTPC, dcaITSTPCTmp[0], dcaITSTPCTmp[1]); - } else if (track.hasASideClustersOnly()) { + } else if (trackFull.hasASideClustersOnly()) { mBufferVals[iThread].back().emplace_back_ITSTPC(Side::A, tglBin, phiBin, qPtBin, multBin, dca[0], dcaZFromDeltaTime, dcarW, dedxRatioqTot, dedxRatioqMax, sqrtChi2TPC, nClTPC, dcaITSTPCTmp[0], dcaITSTPCTmp[1]); } } @@ -1168,7 +1192,7 @@ class TPCTimeSeries : public Task writeData = o2::math_utils::Tsallis::downsampleTsallisCharged(tracksTPC[iTrk].getPt(), factorPt, mSqrt, weight, distr(mGenerator)); } if (writeData) { - auto clusterMask = makeClusterBitMask(track); + auto clusterMask = makeClusterBitMask(trackFull); const auto& trkOrig = tracksTPC[iTrk]; const bool isNearestVtx = (idxITSTPC.back() == -1); // is nearest vertex in case no vertex was found const float mx_ITS = hasITSTPC ? tracksITSTPC[idxITSTPC.front()].getX() : -1; @@ -1177,12 +1201,30 @@ class TPCTimeSeries : public Task const int nClITS = idxITSCheck ? tracksITS[idxITSTrack].getNClusters() : -1; const int chi2ITS = idxITSCheck ? tracksITS[idxITSTrack].getChi2() : -1; int typeSide = 2; // A- and C-Side cluster - if (track.hasASideClustersOnly()) { + if (trackFull.hasASideClustersOnly()) { typeSide = 0; - } else if (track.hasCSideClustersOnly()) { + } else if (trackFull.hasCSideClustersOnly()) { typeSide = 1; } + // check for TOF and propagate TPC track to TOF cluster + bool hasTOFCluster = (std::get<0>(idxTPCTrackToTOFCluster[iTrk]) != -1); + auto tofCl = hasTOFCluster ? tofClusters[std::get<0>(idxTPCTrackToTOFCluster[iTrk])] : o2::tof::Cluster(); + + float tpcYDeltaAtTOF = -999; + float tpcZDeltaAtTOF = -999; + if (hasTOFCluster) { + o2::track::TrackPar trackTmpOut(tracksTPC[iTrk].getParamOut()); + if (trackTmpOut.rotate(o2::math_utils::sector2Angle(tofCl.getSector())) && propagator->propagateTo(trackTmpOut, tofCl.getX(), false, mMaxSnp, mFineStep, mMatType)) { + tpcYDeltaAtTOF = trackTmpOut.getY() - tofCl.getY(); + tpcZDeltaAtTOF = signSide * (o2::tpc::ParameterElectronics::Instance().ZbinWidth * trackFull.getTime0() - vertex.getTimeStamp().getTimeStamp()) * mVDrift - trackTmpOut.getZ() + tofCl.getZ(); + } + } + + // get delta parameter between inner and outer + float deltaTPCParamInOutTgl = trackFull.getTgl() - trackFull.getParamOut().getTgl(); + float deltaTPCParamInOutQPt = trackFull.getQ2Pt() - trackFull.getParamOut().getQ2Pt(); + *mStreamer[iThread] << "treeTimeSeries" // DCAs << "factorPt=" << factorPt @@ -1231,6 +1273,14 @@ class TPCTimeSeries : public Task << "mVDrift=" << mVDrift << "its_flag=" << int(gID) << "sqrtChi2Match=" << chi2Match + // TOF cluster + << "tpcYDeltaAtTOF=" << tpcYDeltaAtTOF + << "tpcZDeltaAtTOF=" << tpcZDeltaAtTOF + << "mDXatTOF=" << std::get<1>(idxTPCTrackToTOFCluster[iTrk]) + << "mDZatTOF=" << std::get<2>(idxTPCTrackToTOFCluster[iTrk]) + // TPC delta param + << "deltaTPCParamInOutTgl=" << deltaTPCParamInOutTgl + << "deltaTPCParamInOutQPt=" << deltaTPCParamInOutQPt << "\n"; } } @@ -1591,7 +1641,7 @@ o2::framework::DataProcessorSpec getTPCTimeSeriesSpec(const bool disableWriter, {"max-dedx-ratio", VariantType::Float, 0.3f, {"Maximum absolute log(dedx(pion)/dedx) ratio"}}, {"max-dedx-region-ratio", VariantType::Float, 0.5f, {"Maximum absolute log(dedx(region)/dedx) ratio"}}, {"sample-unbinned-tsallis", VariantType::Bool, false, {"Perform sampling of unbinned data based on Tsallis function"}}, - {"sampling-factor", VariantType::Float, 0.1f, {"Sampling factor in case sample-unbinned-tsallis is used"}}, + {"sampling-factor", VariantType::Float, 0.001f, {"Sampling factor in case sample-unbinned-tsallis is used"}}, {"out-file-unbinned", VariantType::String, "time_series_tracks.root", {"name of the output file for the unbinned data"}}}}; } diff --git a/Detectors/TPC/workflow/src/TPCTimeSeriesWriterSpec.cxx b/Detectors/TPC/workflow/src/TPCTimeSeriesWriterSpec.cxx index 0160f4a760969..c113a537cef83 100644 --- a/Detectors/TPC/workflow/src/TPCTimeSeriesWriterSpec.cxx +++ b/Detectors/TPC/workflow/src/TPCTimeSeriesWriterSpec.cxx @@ -33,8 +33,8 @@ DataProcessorSpec getTPCTimeSeriesWriterSpec() return MakeRootTreeWriterSpec("tpc-time-series-writer", "o2_timeseries_tpc.root", "treeTimeSeries", - BranchDefinition{InputSpec{"timeseries", o2::header::gDataOriginTPC, getDataDescriptionTimeSeries(), 0}, "TimeSeries", 1}, - BranchDefinition{InputSpec{"itpctfid", o2::header::gDataOriginTPC, getDataDescriptionTPCTimeSeriesTFId(), 0}, "tfID", 1})(); + BranchDefinition{InputSpec{"timeseries", o2::header::gDataOriginTPC, getDataDescriptionTimeSeries(), 0, Lifetime::Sporadic}, "TimeSeries", 1}, + BranchDefinition{InputSpec{"itpctfid", o2::header::gDataOriginTPC, getDataDescriptionTPCTimeSeriesTFId(), 0, Lifetime::Sporadic}, "tfID", 1})(); } } // namespace tpc diff --git a/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx b/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx index 8ea4d67572976..4054eb97158c8 100644 --- a/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx +++ b/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx @@ -586,16 +586,7 @@ bool TRDGlobalTracking::refitITSTPCTRDTrack(TrackTRD& trk, float timeTRD, o2::gl return false; } auto posEnd = trk.getXYZGlo(); - // account path integrals - float dX = posEnd.x() - posStart.x(), dY = posEnd.y() - posStart.y(), dZ = posEnd.z() - posStart.z(), d2XY = dX * dX + dY * dY; - if (std::abs(o2::base::Propagator::Instance()->getNominalBz()) > 0.01) { // circular arc = 2*R*asin(dXY/2R) - float b[3]; - o2::math_utils::Point3D posAv(0.5 * (posEnd.x() + posStart.x()), 0.5 * (posEnd.y() + posStart.y()), 0.5 * (posEnd.z() + posStart.z())); - propagator->getFieldXYZ(posAv, b); - float curvH = std::abs(0.5f * trk.getCurvature(b[2])), arcXY = 1. / curvH * std::asin(curvH * std::sqrt(d2XY)); - d2XY = arcXY * arcXY; - } - auto lInt = std::sqrt(d2XY + dZ * dZ); + auto lInt = propagator->estimateLTIncrement(trk, posStart, posEnd); trk.getLTIntegralOut().addStep(lInt, trk.getP2Inv()); // trk.getLTIntegralOut().addX2X0(lInt * mTPCmeanX0Inv); // do we need to account for the material budget here? probably @@ -685,16 +676,7 @@ bool TRDGlobalTracking::refitTPCTRDTrack(TrackTRD& trk, float timeTRD, o2::globa } auto posEnd = trk.getXYZGlo(); - // account path integrals - float dX = posEnd.x() - posStart.x(), dY = posEnd.y() - posStart.y(), dZ = posEnd.z() - posStart.z(), d2XY = dX * dX + dY * dY; - if (std::abs(o2::base::Propagator::Instance()->getNominalBz()) > 0.01) { // circular arc = 2*R*asin(dXY/2R) - float b[3]; - o2::math_utils::Point3D posAv(0.5 * (posEnd.x() + posStart.x()), 0.5 * (posEnd.y() + posStart.y()), 0.5 * (posEnd.z() + posStart.z())); - propagator->getFieldXYZ(posAv, b); - float curvH = std::abs(0.5f * trk.getCurvature(b[2])), arcXY = 1. / curvH * std::asin(curvH * std::sqrt(d2XY)); - d2XY = arcXY * arcXY; - } - auto lInt = std::sqrt(d2XY + dZ * dZ); + auto lInt = propagator->estimateLTIncrement(trk, posStart, posEnd); trk.getLTIntegralOut().addStep(lInt, trk.getP2Inv()); // trk.getLTIntegralOut().addX2X0(lInt * mTPCmeanX0Inv); // do we need to account for the material budget here? probably? diff --git a/Detectors/Upgrades/ITS3/base/CMakeLists.txt b/Detectors/Upgrades/ITS3/base/CMakeLists.txt index eb25cfe651ea1..1843279216345 100644 --- a/Detectors/Upgrades/ITS3/base/CMakeLists.txt +++ b/Detectors/Upgrades/ITS3/base/CMakeLists.txt @@ -11,13 +11,8 @@ o2_add_library(ITS3Base SOURCES src/MisalignmentParameter.cxx - src/SegmentationSuperAlpide.cxx - src/SuperAlpideParams.cxx - src/DescriptorInnerBarrelITS3Param.cxx PUBLIC_LINK_LIBRARIES O2::CommonConstants O2::MathUtils O2::DetectorsBase) o2_target_root_dictionary(ITS3Base HEADERS include/ITS3Base/SegmentationSuperAlpide.h - include/ITS3Base/MisalignmentParameter.h - include/ITS3Base/SuperAlpideParams.h - include/ITS3Base/DescriptorInnerBarrelITS3Param.h) + include/ITS3Base/MisalignmentParameter.h) diff --git a/Detectors/Upgrades/ITS3/base/include/ITS3Base/DescriptorInnerBarrel.h b/Detectors/Upgrades/ITS3/base/include/ITS3Base/DescriptorInnerBarrel.h new file mode 100644 index 0000000000000..e69de29bb2d1d diff --git a/Detectors/Upgrades/ITS3/base/include/ITS3Base/DescriptorInnerBarrelITS3Param.h b/Detectors/Upgrades/ITS3/base/include/ITS3Base/DescriptorInnerBarrelITS3Param.h deleted file mode 100644 index b6a44d020c414..0000000000000 --- a/Detectors/Upgrades/ITS3/base/include/ITS3Base/DescriptorInnerBarrelITS3Param.h +++ /dev/null @@ -1,55 +0,0 @@ -// 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. - -/// \file DescriptorInnerBarrelITS3Param.h -/// \brief Definition of the DescriptorInnerBarrelITS3Param class - -#ifndef ALICEO2_ITS3_DESCRIPTORINNERBARRELITS3PARAM_H -#define ALICEO2_ITS3_DESCRIPTORINNERBARRELITS3PARAM_H - -#include "CommonUtils/ConfigurableParam.h" -#include "CommonUtils/ConfigurableParamHelper.h" - -namespace o2 -{ -namespace its3 -{ -/** - ** a parameter class/struct to keep the settings of - ** the DescriptorInnerBarrelITS3 class - ** allow the user to modify them - **/ - -enum class ITS3Version { - None = 0, /* none */ - ThreeLayersNoDeadZones = 1, /* three layers without dead zones */ - ThreeLayers = 2, /* three layers with dead zones */ - FourLayers = 3, /* four layers with dead zones */ - ThreeLayersDeadZonesFirstOnly = 4 /* three layers with dead zones only in first layer */ -}; - -struct DescriptorInnerBarrelITS3Param : public o2::conf::ConfigurableParamHelper { - ITS3Version mVersion = ITS3Version::None; - int mBuildLevel{0}; - double mGapY[4] = {0.f, 0.f, 0.f, 0.f}; - double mGapPhi[4] = {0.1f, 0.1f, 0.1f, 0.1f}; - double mRadii[4] = {1.8f, 2.4f, 3.0f, 6.0f}; - double mLength{26.f}; - double mGapXDirection4thLayer{0.f}; - double mAddMaterial3rdLayer{0.f}; - std::string const& getITS3LayerConfigString() const; - O2ParamDef(DescriptorInnerBarrelITS3Param, "DescriptorInnerBarrelITS3"); -}; - -} // namespace its3 -} // end namespace o2 - -#endif // ALICEO2_ITS3_DESCRIPTORINNERBARRELITS3PARAM_H \ No newline at end of file diff --git a/Detectors/Upgrades/ITS3/base/include/ITS3Base/GeometryTGeo.h b/Detectors/Upgrades/ITS3/base/include/ITS3Base/GeometryTGeo.h new file mode 100644 index 0000000000000..89331597b89fc --- /dev/null +++ b/Detectors/Upgrades/ITS3/base/include/ITS3Base/GeometryTGeo.h @@ -0,0 +1,69 @@ +// 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. + +/// \file GeometryTGeo.h +/// \brief Definition of the GeometryTGeo class for ITS3 +/// \author felix.schlepper@cern.ch + +#ifndef ALICEO2_ITS3_GEOMETRYTGEO_H_ +#define ALICEO2_ITS3_GEOMETRYTGEO_H_ + +#include +#include +#include +#include "DetectorsBase/GeometryManager.h" +#include "DetectorsCommonDataFormats/DetID.h" +#include "DetectorsCommonDataFormats/DetMatrixCache.h" +#include "MathUtils/Utils.h" + +namespace o2::its3 +{ +class GeometryTGeo : public o2::detectors::DetMatrixCache +{ + public: + using Mat3D = o2::math_utils::Transform3D; + using DetMatrixCache::getMatrixL2G; + using DetMatrixCache::getMatrixT2G; + using DetMatrixCache::getMatrixT2GRot; + using DetMatrixCache::getMatrixT2L; + using o2::detectors::DetMatrixCache::fillMatrixCache; + + GeometryTGeo() = default; + GeometryTGeo(const o2::detectors::DetID& detid) : o2::detectors::DetMatrixCache(detid) {} + + GeometryTGeo(const GeometryTGeo& src) = delete; + GeometryTGeo(GeometryTGeo&& src) = delete; + GeometryTGeo& operator=(const GeometryTGeo& geom) = delete; + ~GeometryTGeo() final = default; + + static GeometryTGeo* Instance() + { + // get (create if needed) a unique instance of the object +#ifdef GPUCA_STANDALONE + return nullptr; // TODO: DR: Obviously wrong, but to make it compile for now +#else + if (!mInstance) { + mInstance = std::unique_ptr(new GeometryTGeo()); + } + return mInstance.get(); +#endif + } + + private: +#ifndef GPUCA_STANDALONE + static std::unique_ptr mInstance; ///< singletone instance +#endif + + ClassDefOverride(GeometryTGeo, 0); +}; +} // namespace o2::its3 + +#endif diff --git a/Detectors/Upgrades/ITS3/base/include/ITS3Base/SegmentationSuperAlpide.h b/Detectors/Upgrades/ITS3/base/include/ITS3Base/SegmentationSuperAlpide.h index fb888feb359ea..d178e88e66de4 100644 --- a/Detectors/Upgrades/ITS3/base/include/ITS3Base/SegmentationSuperAlpide.h +++ b/Detectors/Upgrades/ITS3/base/include/ITS3Base/SegmentationSuperAlpide.h @@ -11,16 +11,15 @@ /// \file SegmentationSuperAlpide.h /// \brief Definition of the SegmentationSuperAlpide class +/// \author Fabrizio Grosa +/// \author felix.schlepper@cern.ch #ifndef ALICEO2_ITS3_SEGMENTATIONSUPERALPIDE_H_ #define ALICEO2_ITS3_SEGMENTATIONSUPERALPIDE_H_ -#include #include "MathUtils/Cartesian.h" #include "CommonConstants/MathConstants.h" -#include "ITS3Base/SuperAlpideParams.h" -#include "ITS3Base/DescriptorInnerBarrelITS3Param.h" -#include +#include "ITSBase/Specs.h" namespace o2 { @@ -28,40 +27,37 @@ namespace its3 { /// Segmentation and response for pixels in ITS3 upgrade +template class SegmentationSuperAlpide { + // This class defines the segmenation of the pixelArray in the tile. We define + // two coordinate systems, one width x,z detector local coordianates (cm) and + // the more natural row,col layout: Also all the transformation between these + // two. The class provides the transformation from the tile to TGeo + // coordinates. + + // row,col=0 + // | + // v + // x----------------------x + // | | | + // | | | + // | | | + // | | | ^ x + // | | | | + // | | | | + // | | | | + // |-----------X----------| X marks (x,z)=(0,0) X----> z + // | | | + // | | | + // | | | + // | | | + // | | | + // | | | + // x----------------------x + public: - SegmentationSuperAlpide(int layer, float pitchCol, float pitchRow, float detThick, double length, const double radii[4]) : mLayer{layer}, mPitchCol{pitchCol}, mPitchRow{pitchRow}, mDetectorLayerThickness{detThick}, mLength{length} - { - for (int iL{0}; iL < 4; ++iL) { - mRadii[iL] = radii[iL]; - } - print(); - } - SegmentationSuperAlpide(int layer = 0) : SegmentationSuperAlpide(layer, SuperAlpideParams::Instance().mPitchCol, SuperAlpideParams::Instance().mPitchRow, SuperAlpideParams::Instance().mDetectorThickness, DescriptorInnerBarrelITS3Param::Instance().mLength, DescriptorInnerBarrelITS3Param::Instance().mRadii) {} - - double mRadii[4] = {1.8f, 2.4f, 3.0f, 6.0f}; ///< radii for different layers - const double mLength; ///< chip length - const int mLayer; ///< chip layer - const float mPitchCol; ///< pixel column size at the external surface - const float mPitchRow; ///< pixel row size at the external surface - const float mDetectorLayerThickness; ///< detector thickness - const int mNCols{static_cast(std::floor(mLength / mPitchCol / 2) * 2)}; ///< number of columns (we force to have an even number) - const int mNRows{static_cast(std::floor(double(mRadii[mLayer] + mDetectorLayerThickness - mSensorLayerThicknessEff / 2.) * double(constants::math::PI) / double(mPitchRow) / 2 * (mLayer == 3 ? 0.5 : 1.)) * 2)}; ///< number of rows (we force to have an even number) - const int mNPixels{mNRows * mNCols}; ///< total number of pixels - static constexpr float mPassiveEdgeReadOut = 0.; ///< width of the readout edge (Passive bottom) - static constexpr float mPassiveEdgeTop = 0.; ///< Passive area on top - static constexpr float mPassiveEdgeSide = 0.; ///< width of Passive area on left/right of the sensor - const float mActiveMatrixSizeCols{mPitchCol * mNCols}; ///< Active size along columns - const float mActiveMatrixSizeRows{mPitchRow * mNRows}; ///< Active size along rows - - // effective thickness of sensitive layer, accounting for charge collection non-uniformity, https://alice.its.cern.ch/jira/browse/AOC-46 - static constexpr float mSensorLayerThicknessEff = 28.e-4; ///< effective thickness of sensitive part - static constexpr float mSensorLayerThickness = 30.e-4; ///< physical thickness of sensitive part - const float mSensorSizeCols{mActiveMatrixSizeCols + mPassiveEdgeSide + mPassiveEdgeSide}; ///< SensorSize along columns - const float mSensorSizeRows{mActiveMatrixSizeRows + mPassiveEdgeTop + mPassiveEdgeReadOut}; ///< SensorSize along rows - - ~SegmentationSuperAlpide() = default; + SegmentationSuperAlpide(int layer = 0) : mLayer{layer} {} /// Transformation from the curved surface to a flat surface /// It works only if the detector is not rototraslated @@ -73,7 +69,13 @@ class SegmentationSuperAlpide /// the center of the sensitive volume. /// \param yFlat Detector local flat coordinate y in cm with respect to /// the center of the sensitive volume. - void curvedToFlat(float xCurved, float yCurved, float& xFlat, float& yFlat); + void curvedToFlat(value_t xCurved, value_t yCurved, value_t& xFlat, value_t& yFlat) + { + value_t dist = std::hypot(xCurved, yCurved); + yFlat = dist - mEffRadius; + value_t phi = (vaule_t)constants::math::PI / 2 - std::atan2((double)yCurved, (double)xCurved); + xFlat = mEffRadius * phi; + } /// Transformation from the flat surface to a curved surface /// It works only if the detector is not rototraslated @@ -85,145 +87,122 @@ class SegmentationSuperAlpide /// the center of the sensitive volume. /// \param yCurved Detector local curved coordinate y in cm with respect to /// the center of the sensitive volume. - void flatToCurved(float xFlat, float yFlat, float& xCurved, float& yCurved); + void flatToCurved(value_t xFlat, value_t yFlat, value_t& xCurved, value_t& yCurved) + { + value_t dist = yFlat + mEffRadius; + value_t phi = xFlat / dist; + value_t tang = std::tan((value_t)constants::math::PI / 2 - (value_t)phi); + xCurved = (xFlat > 0 ? 1.f : -1.f) * dist / std::sqrt(1 + tang * tang); + yCurved = xCurved * tang; + } /// Transformation from Geant detector centered local coordinates (cm) to /// Pixel cell numbers iRow and iCol. - /// Returns kTRUE if point x,z is inside sensitive volume, kFALSE otherwise. + /// Returns true if point x,z is inside sensitive volume, false otherwise. /// A value of -1 for iRow or iCol indicates that this point is outside of the /// detector segmentation as defined. /// \param float x Detector local coordinate x in cm with respect to /// the center of the sensitive volume. /// \param float z Detector local coordinate z in cm with respect to /// the center of the sensitive volume. - /// \param int iRow Detector x cell coordinate. Has the range 0 <= iRow < mNumberOfRows - /// \param int iCol Detector z cell coordinate. Has the range 0 <= iCol < mNumberOfColumns - bool localToDetector(float x, float z, int& iRow, int& iCol); - /// same but w/o check for row/column range - void localToDetectorUnchecked(float xRow, float zCol, int& iRow, int& iCol); + /// \param int iRow Detector x cell coordinate. + /// \param int iCol Detector z cell coordinate. + bool localToDetector(value_t const xRow, value_t const zCol, int& iRow, int& iCol) const noexcept + { + localToDetectorUnchecked(xRow, zCol, iRow, iCol); + if (!isValid(iRow, iCol)) { + iRow = iCol = -1; + return false; + } + return true; + } - /// Transformation from Detector cell coordiantes to Geant detector centered + // Same as localToDetector w.o. checks. + void localToDetectorUnchecked(value_t const xRow, value_t const zCol, int& iRow, int& iCol) const noexcept + { + namespace cp = constants::pixelarray; + value_t x = cp::length / 2. - xRow; // transformation to upper edge of pixelarray + value_t z = zCol + cp::width / 2.; // transformation to left edge of pixelarray + iRow = std::floor(x / cp::pixel::pitchRow); + iCol = std::floor(z / cp::pixel::pitchCol); + } + + /// Transformation from Detector cell coordinates to Geant detector centered /// local coordinates (cm) - /// \param int iRow Detector x cell coordinate. Has the range 0 <= iRow < mNumberOfRows - /// \param int iCol Detector z cell coordinate. Has the range 0 <= iCol < mNumberOfColumns + /// \param int iRow Detector x cell coordinate. + /// \param int iCol Detector z cell coordinate. /// \param float x Detector local coordinate x in cm with respect to the /// center of the sensitive volume. /// \param float z Detector local coordinate z in cm with respect to the /// center of the sensitive volume. /// If iRow and or iCol is outside of the segmentation range a value of -0.5*Dx() /// or -0.5*Dz() is returned. - bool detectorToLocal(int iRow, int iCol, float& xRow, float& zCol); - bool detectorToLocal(float row, float col, float& xRow, float& zCol); - bool detectorToLocal(float row, float col, math_utils::Point3D& loc); - - // same but w/o check for row/col range - void detectorToLocalUnchecked(int iRow, int iCol, float& xRow, float& zCol); - void detectorToLocalUnchecked(float row, float col, float& xRow, float& zCol); - void detectorToLocalUnchecked(float row, float col, math_utils::Point3D& loc); - - float getFirstRowCoordinate() + bool detectorToLocal(int const iRow, int const iCol, value_t& xRow, value_t& zCol) const noexcept { - return 0.5 * ((mActiveMatrixSizeRows - mPassiveEdgeTop + mPassiveEdgeReadOut) - mPitchRow); + detectorToLocalUnchecked(iRow, iCol, xRow, zCol); + if (!isValid(xRow, zCol)) { + return false; + } + return true; } - float getFirstColCoordinate() + + // Same as detectorToLocal w.o. checks. + // We position ourself in the middle of the pixel. + void detectorToLocalUnchecked(int iRow, int iCol, value_t& xRow, value_t& zCol) const noexcept { - return 0.5 * (mPitchCol - mActiveMatrixSizeCols); + namespace cp = constants::pixelarray; + xRow = -(iRow - 0.5) * cp::pixel::pitchRow + cp::length / 2.; + zCol = -(iCol + 0.5) * cp::pixel::pitchCol - cp::width / 2.; } - void print(); - - ClassDefNV(SegmentationSuperAlpide, 2); // Segmentation class upgrade pixels -}; // namespace its3 - -inline void SegmentationSuperAlpide::curvedToFlat(float xCurved, float yCurved, float& xFlat, float& yFlat) -{ - float dist = std::sqrt(xCurved * xCurved + yCurved * yCurved); - yFlat = dist - (mRadii[mLayer] + mDetectorLayerThickness - mSensorLayerThickness / 2.); - float phi = (double)constants::math::PI / 2 - std::atan2((double)yCurved, (double)xCurved); - xFlat = (mRadii[mLayer] + mDetectorLayerThickness - mSensorLayerThickness / 2.) * phi; // we bring everything to the upper surface to avoid effects due to the different length of upper and lower surfaces -} - -inline void SegmentationSuperAlpide::flatToCurved(float xFlat, float yFlat, float& xCurved, float& yCurved) -{ - float dist = yFlat + mRadii[mLayer] + mDetectorLayerThickness - mSensorLayerThickness / 2.; - float phi = xFlat / dist; - float tang = std::tan((double)constants::math::PI / 2 - (double)phi); - xCurved = (xFlat > 0 ? 1.f : -1.f) * dist / std::sqrt(1 + tang * tang); - yCurved = xCurved * tang; -} - -inline void SegmentationSuperAlpide::localToDetectorUnchecked(float xRow, float zCol, int& iRow, int& iCol) -{ - // convert to row/col w/o over/underflow check - xRow = 0.5 * (mActiveMatrixSizeRows - mPassiveEdgeTop + mPassiveEdgeReadOut) - xRow; // coordinate wrt top edge of Active matrix - zCol += 0.5 * mActiveMatrixSizeCols; // coordinate wrt left edge of Active matrix - iRow = std::floor(xRow / mPitchRow); - iCol = std::floor(zCol / mPitchCol); - if (xRow < 0) { - iRow -= 1; - } - if (zCol < 0) { - iCol -= 1; + bool detectorToLocal(float row, float col, float& xRow, float& zCol) + { + return detectorToLocal(static_cast(row), static_cast(col), xRow, zCol); } -} -inline bool SegmentationSuperAlpide::localToDetector(float xRow, float zCol, int& iRow, int& iCol) -{ - // convert to row/col - xRow = 0.5 * (mActiveMatrixSizeRows - mPassiveEdgeTop + mPassiveEdgeReadOut) - xRow; // coordinate wrt left edge of Active matrix - zCol += 0.5 * mActiveMatrixSizeCols; // coordinate wrt bottom edge of Active matrix - if (xRow < 0 || xRow >= mActiveMatrixSizeRows || zCol < 0 || zCol >= mActiveMatrixSizeCols) { - iRow = iCol = -1; - return false; + void detectorToLocalUnchecked(float row, float col, float& xRow, float& zCol) + { + detectorToLocalUnchecked(static_cast(row), static_cast(col), xRow, zCol); } - iRow = std::floor(xRow / mPitchRow); - iCol = std::floor(zCol / mPitchCol); - return true; -} - -inline void SegmentationSuperAlpide::detectorToLocalUnchecked(int iRow, int iCol, float& xRow, float& zCol) -{ - xRow = getFirstRowCoordinate() - iRow * mPitchRow; - zCol = iCol * mPitchCol + getFirstColCoordinate(); -} - -inline void SegmentationSuperAlpide::detectorToLocalUnchecked(float row, float col, float& xRow, float& zCol) -{ - xRow = getFirstRowCoordinate() - row * mPitchRow; - zCol = col * mPitchCol + getFirstColCoordinate(); -} -inline void SegmentationSuperAlpide::detectorToLocalUnchecked(float row, float col, math_utils::Point3D& loc) -{ - loc.SetCoordinates(getFirstRowCoordinate() - row * mPitchRow, 0.f, col * mPitchCol + getFirstColCoordinate()); -} - -inline bool SegmentationSuperAlpide::detectorToLocal(int iRow, int iCol, float& xRow, float& zCol) -{ - if (iRow < 0 || iRow >= mNRows || iCol < 0 || iCol >= mNCols) { - return false; + bool detectorToLocal(float row, float col, math_utils::Point3D& loc) + { + float xRow, zCol; + if (detectorToLocal(row, col, xRow, zCol)) { + return false; + } + loc.SetCoordinates(xRow, 0., zCol); + return true; + } + void detectorToLocalUnchecked(float row, float col, math_utils::Point3D& loc) + { + float xRow, zCol; + detectorToLocalUnchecked(row, col, xRow, zCol); + loc.SetCoordinates(xRow, 0., zCol); } - detectorToLocalUnchecked(iRow, iCol, xRow, zCol); - return true; -} -inline bool SegmentationSuperAlpide::detectorToLocal(float row, float col, float& xRow, float& zCol) -{ - if (row < 0 || row >= mNRows || col < 0 || col >= mNCols) { - return false; + private: + bool isValid(value_t const xRow, value_t const zCol) const noexcept + { + namespace cp = constants::pixelarray; + if (xRow < 0. || xRow >= cp::length || zCol < 0. || zCol >= cp::width) { + return false; + } + return true; } - detectorToLocalUnchecked(row, col, xRow, zCol); - return true; -} -inline bool SegmentationSuperAlpide::detectorToLocal(float row, float col, math_utils::Point3D& loc) -{ - if (row < 0 || row >= mNRows || col < 0 || col >= mNCols) { - return false; + bool isValid(int const iRow, int const iCol) const noexcept + { + namespace cp = constants::pixelarray; + if (iRow < 0 || iRow >= cp::nRows || iCol < 0 || iCol >= cp::nCols) { + return false; + } + return true; } - detectorToLocalUnchecked(row, col, loc); - return true; -} + + const int mLayer; ///< chip layer + const value_t mEffRadius{constants::radii[mLayer] + constants::thickness / 2.}; +}; } // namespace its3 } // namespace o2 diff --git a/Detectors/Upgrades/ITS3/base/include/ITS3Base/Specs.h b/Detectors/Upgrades/ITS3/base/include/ITS3Base/Specs.h new file mode 100644 index 0000000000000..4ef824f73b179 --- /dev/null +++ b/Detectors/Upgrades/ITS3/base/include/ITS3Base/Specs.h @@ -0,0 +1,101 @@ +// 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. + +/// \file Specs.h +/// \brief TDR specs of ITS3 +/// \author felix.schlepper@cern.ch + +#include "Rtypes.h" + +namespace o2::its3 +{ +namespace constants +{ +constexpr double cm{1e+2}; // This is the default unit of TGeo so we use this as scale +constexpr double mu{1e-6 * cm}; +constexpr double mm{1e-3 * cm}; +namespace pixelarray +{ +constexpr double length{9.197 * mm}; +constexpr double width{3.571 * mm}; +constexpr EColor color{kGreen}; +constexpr unsigned int nCols{156}; +constexpr unsigned int nRows{440}; +constexpr unsigned int nPixels{nCols * nRows}; +namespace pixel +{ +constexpr double pitchCol{width / static_cast(nCols)}; +constexpr double pitchRow{length / static_cast(nRows)}; +} // namespace pixel +} // namespace pixelarray +namespace tile +{ +namespace biasing +{ +constexpr double length{0.06 * mm}; +constexpr double width{3.571 * mm}; +constexpr EColor color{kYellow}; +static_assert(width == pixelarray::width); +} // namespace biasing +namespace powerswitches +{ +constexpr double length{9.257 * mm}; +constexpr double width{0.02 * mm}; +constexpr double z{pixelarray::width}; +constexpr EColor color{kBlue}; +} // namespace powerswitches +namespace readout +{ +constexpr double length{0.525 * mm}; +constexpr double width{3.591 * mm}; +constexpr EColor color{kMagenta}; +static_assert(width == (biasing::width + powerswitches::width)); +} // namespace readout +constexpr double width{readout::width}; +constexpr double length{powerswitches::length + readout::length}; +} // namespace tile +namespace rsu +{ +namespace databackbone +{ +constexpr double length{9.782 * mm}; +constexpr double width{0.06 * mm}; +constexpr EColor color{kRed}; +} // namespace databackbone +constexpr double length{19.564 * mm}; +constexpr double width{21.666 * mm}; +} // namespace rsu +namespace segment +{ +constexpr double length{rsu::length}; +namespace lec +{ +constexpr double length{segment::length}; +constexpr double width{4.5 * mm}; +constexpr EColor color{kCyan}; +} // namespace lec +namespace rec +{ +constexpr double length{segment::length}; +constexpr double width{1.5 * mm}; +constexpr EColor color{kCyan}; +} // namespace rec +constexpr unsigned int nRSUs{12}; +constexpr double width{nRSUs * rsu::width + lec::width + rec::width}; +} // namespace segment +constexpr unsigned int nLayers{3}; +constexpr std::array radii{19 * mm, 25.2 * mm, 31.5 * mm}; // middle radius e.g. inner radius+thickness/2. +constexpr double equatorialGap{1 * mm}; +constexpr std::array nSegments{3, 4, 5}; +constexpr double thickness{50 * mu}; +constexpr double effThickness{66 * mu}; +} // namespace constants +} // namespace o2::its3 diff --git a/Detectors/Upgrades/ITS3/base/include/ITS3Base/SuperAlpideParams.h b/Detectors/Upgrades/ITS3/base/include/ITS3Base/SuperAlpideParams.h deleted file mode 100644 index d3971efcd2431..0000000000000 --- a/Detectors/Upgrades/ITS3/base/include/ITS3Base/SuperAlpideParams.h +++ /dev/null @@ -1,38 +0,0 @@ -// 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. - -// \brief Collect all possible configurable parameters for Super ALPIDE chips - -#ifndef O2_SUPER_ALPIDE_PARAMS_H -#define O2_SUPER_ALPIDE_PARAMS_H - -#include "CommonUtils/ConfigurableParam.h" -#include "CommonUtils/ConfigurableParamHelper.h" - -namespace o2 -{ -namespace its3 -{ - -/// Segmentation parameters for Super ALPIDE chips -struct SuperAlpideParams : public o2::conf::ConfigurableParamHelper { - float mPitchCol = 29.24e-4; ///< Pixel column size (cm) //FIXME: proxy value to get same resolution as ITS2 given incorrect sensor response - float mPitchRow = 26.88e-4; ///< Pixel row size (cm) //FIXME: proxy value to get same resolution as ITS2 given incorrect sensor response - float mDetectorThickness = 50.e-4; ///< Detector thickness (cm) - - // boilerplate - O2ParamDef(SuperAlpideParams, "SuperAlpideParams"); -}; - -} // namespace its3 -} // namespace o2 - -#endif diff --git a/Detectors/Upgrades/ITS3/base/src/DescriptorInnerBarrel.cxx b/Detectors/Upgrades/ITS3/base/src/DescriptorInnerBarrel.cxx new file mode 100644 index 0000000000000..e69de29bb2d1d diff --git a/Detectors/Upgrades/ITS3/base/src/DescriptorInnerBarrelITS3Param.cxx b/Detectors/Upgrades/ITS3/base/src/DescriptorInnerBarrelITS3Param.cxx deleted file mode 100644 index 2e2d8df83a5a2..0000000000000 --- a/Detectors/Upgrades/ITS3/base/src/DescriptorInnerBarrelITS3Param.cxx +++ /dev/null @@ -1,34 +0,0 @@ -// 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. - -/// \file DescriptorInnerBarrelITS3Param.h -/// \brief Implementation of the DescriptorInnerBarrelITS3Param class - -#include "ITS3Base/DescriptorInnerBarrelITS3Param.h" -O2ParamImpl(o2::its3::DescriptorInnerBarrelITS3Param); - -namespace o2 -{ -namespace its3 -{ - -namespace -{ -static const std::string confstrings[5] = {"", "ThreeLayersNoDeadZones", "ThreeLayers", "FourLayers", "FiveLayers"}; -} - -std::string const& DescriptorInnerBarrelITS3Param::getITS3LayerConfigString() const -{ - return confstrings[(int)mVersion]; -} - -} // namespace its3 -} // end namespace o2 diff --git a/Detectors/Upgrades/ITS3/base/src/SuperAlpideParams.cxx b/Detectors/Upgrades/ITS3/base/src/GeometryTGeo.cxx similarity index 66% rename from Detectors/Upgrades/ITS3/base/src/SuperAlpideParams.cxx rename to Detectors/Upgrades/ITS3/base/src/GeometryTGeo.cxx index cbaa8b06cc15f..b530b03a5cdf1 100644 --- a/Detectors/Upgrades/ITS3/base/src/SuperAlpideParams.cxx +++ b/Detectors/Upgrades/ITS3/base/src/GeometryTGeo.cxx @@ -9,5 +9,14 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -#include "ITS3Base/SuperAlpideParams.h" -O2ParamImpl(o2::its3::SuperAlpideParams); +/// \file GeometryTGeo.cxx +/// \brief Implementation of the GeometryTGeo class +/// \author felix.schlepper@cern.ch + +#include "ITS3Base/GeometryTGeo.h" + +#include "Rtypes.h" + +ClassImp(o2::its::GeometryTGeo); + +std::unique_ptr GeometryTGeo::mInstance; diff --git a/Detectors/Upgrades/ITS3/base/src/ITS3BaseLinkDef.h b/Detectors/Upgrades/ITS3/base/src/ITS3BaseLinkDef.h index 854b296435b69..89e8790bb58ce 100644 --- a/Detectors/Upgrades/ITS3/base/src/ITS3BaseLinkDef.h +++ b/Detectors/Upgrades/ITS3/base/src/ITS3BaseLinkDef.h @@ -15,11 +15,6 @@ #pragma link off all classes; #pragma link off all functions; -#pragma link C++ class o2::its3::SegmentationSuperAlpide + ; #pragma link C++ class o2::its3::MisalignmentParameter + ; -#pragma link C++ class o2::its3::SuperAlpideParams + ; -#pragma link C++ class o2::its3::DescriptorInnerBarrelITS3Param + ; -#pragma link C++ class o2::conf::ConfigurableParamHelper < o2::its3::SuperAlpideParams> + ; -#pragma link C++ class o2::conf::ConfigurableParamHelper < o2::its3::DescriptorInnerBarrelITS3Param> + ; #endif diff --git a/Detectors/Upgrades/ITS3/base/src/SegmentationSuperAlpide.cxx b/Detectors/Upgrades/ITS3/base/src/SegmentationSuperAlpide.cxx deleted file mode 100644 index 8245d4c205d14..0000000000000 --- a/Detectors/Upgrades/ITS3/base/src/SegmentationSuperAlpide.cxx +++ /dev/null @@ -1,26 +0,0 @@ -// 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. - -/// \file SegmentationSuperAlpide.cxx -/// \brief Implementation of the SegmentationSuperAlpide class - -#include "ITS3Base/SegmentationSuperAlpide.h" -#include - -ClassImp(o2::its3::SegmentationSuperAlpide); - -void o2::its3::SegmentationSuperAlpide::print() -{ - LOGP(debug, "SegmentationSuperAlpide:"); - LOGP(debug, "Layer {}: Active/Total size {:.2f}/{:.2f} (rows) x {:.2f}/{:.2f}", mLayer, mActiveMatrixSizeRows, mSensorSizeRows, mActiveMatrixSizeCols, mSensorSizeCols); - LOGP(debug, "Pixel size: {:.2f} (along {} rows) {:.2f} (along {} columns) microns", mPitchRow * 1e4, mNRows, mPitchCol * 1e4, mNCols); - LOGP(debug, "Passive edges: bottom: {:.6f}, top: {:.6f}, left/right: {:.6f} microns", mPassiveEdgeReadOut * 1e4, mPassiveEdgeTop * 1e4, mPassiveEdgeSide * 1e4); -} diff --git a/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/Detector.h b/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/Detector.h new file mode 100644 index 0000000000000..e69de29bb2d1d diff --git a/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/ITS3Layer.h b/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/ITS3Layer.h index c97ff93a01e72..1e0465b599d57 100644 --- a/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/ITS3Layer.h +++ b/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/ITS3Layer.h @@ -12,81 +12,62 @@ /// \file ITS3Layer.h /// \brief Definition of the ITS3Layer class /// \author Fabrizio Grosa +/// \author felix.schlepper@cern.ch #ifndef ALICEO2_ITS3_ITS3LAYER_H #define ALICEO2_ITS3_ITS3LAYER_H -#include +#include class TGeoVolume; -namespace o2 +namespace o2::its3 { -namespace its3 -{ -/// This class defines the Geometry for the ITS3 using TGeo. This is a work class used -/// to study different configurations during the development of the new ITS structure -class ITS3Layer : public TObject + +/// This class defines the Geometry for the ITS3 using TGeo. +class ITS3Layer { + // The hierarchy will be the following: + // ITS2 -> ITS3 + // --------------------------------- + // Sensor PixelArray + // Chip Tile + // Module RSU + // HalfStave Segment + // Stave Chip + // HalfBarrel CarbonForm + // Layer Layer public: - // Default constructor - ITS3Layer() = default; - - /// Constructor setting layer number - ITS3Layer(int lay); - - /// Copy constructor - ITS3Layer(const ITS3Layer&) = default; - - /// Assignment operator - ITS3Layer& operator=(const ITS3Layer&) = default; + // Create one layer of ITS3 and attach it to the motherVolume. + void createLayer(TGeoVolume* motherVolume, int layer = 0); - /// Default destructor - ~ITS3Layer() override; - - void createLayer(TGeoVolume* motherVolume, double radiusBetweenLayer); - void createLayerWithDeadZones(TGeoVolume* motherVolume, double radiusBetweenLayer); - void create4thLayer(TGeoVolume* motherVolume); - void createCarbonFoamStructure(TGeoVolume* motherVolume, double deltaR, bool fourthLayer = false); - - void setChipThick(double thick) { mChipThickness = thick; } - void setLayerRadius(double radius) { mRadius = radius; } - void setLayerZLen(double zLen) { mZLen = zLen; } - void setGapBetweenEmispheres(double gap) { mGapY = gap; } - void setGapBetweenEmispheresInPhi(double gap) { mGapPhi = gap; } - void setChipID(int chipID) { mChipTypeID = chipID; } - void setFringeChipWidth(double fringeChipWidth) { mFringeChipWidth = fringeChipWidth; } - void setMiddleChipWidth(double middleChipWidth) { mMiddleChipWidth = middleChipWidth; } - void setNumSubSensorsHalfLayer(int numSubSensorsHalfLayer) { mNumSubSensorsHalfLayer = numSubSensorsHalfLayer; } - void setHeightStripFoam(double heightStripFoam) { mHeightStripFoam = heightStripFoam; } - void setLengthSemiCircleFoam(double lengthSemiCircleFoam) { mLengthSemiCircleFoam = lengthSemiCircleFoam; } - void setThickGluedFoam(double thickGluedFoam) { mThickGluedFoam = thickGluedFoam; } - void setGapXDirection(double gapXDirection) { mGapXDirection = gapXDirection; } - void setBuildLevel(int buildLevel) { mBuildLevel = buildLevel; } - void setAdditionalMaterial(double addMat) { mAddMaterial = addMat; } + private: + void init(); + void createPixelArray(); + void createTile(); + void createRSU(); + void createSegment(); + void createChip(); + void createCarbonForm(); + void createLayerImpl(); private: - int mLayerNumber{0}; //! layer number - double mChipThickness{0.}; //! chip thickness - double mSensorThickness{0.}; //! sensor thickness - double mRadius{0.}; //! radius of layer - double mZLen{0.}; //! length of a layer - double mGapY{0.}; //! gap between emispheres in Y direction - double mGapPhi{0.}; //! gap between emispheres in phi (distance in Y direction) - int mChipTypeID{0}; //! chip ID - double mFringeChipWidth{0.}; //! fringe chip width - double mMiddleChipWidth{0.}; //! middle chip width - int mNumSubSensorsHalfLayer{0}; //! num of subsensors in half layer - double mHeightStripFoam{0.}; //! strip foam height - double mLengthSemiCircleFoam{0.}; //! semi-circle foam length - double mThickGluedFoam{0.}; //! glued foam thickness - double mGapXDirection{0.}; //! gap between quarter layer(only for layer 4) - double mAddMaterial{0.}; //! additional material to mimic services - int mBuildLevel{0}; //! build level for material budget studies + uint8_t mNLayer{0}; // Layer number + double mR{0}; // Middle Radius + double mRmin{}; // Minimum Radius + double mRmax{0}; // Maximum Radius + + // Individual Pieces + TGeoVolume* mPixelArray{nullptr}; + TGeoVolumeAssembly* mTile{nullptr}; + TGeoVolumeAssembly* mRSU{nullptr}; + TGeoVolumeAssembly* mSegment{nullptr}; + TGeoVolumeAssembly* mChip{nullptr}; + TGeoVolumeAssembly* mCarbonForm{nullptr}; + TGeoVolumeAssembly* mLayer{nullptr}; - ClassDefOverride(ITS3Layer, 0); // ITS3 geometry -}; // namespace its3 -} // namespace its3 -} // namespace o2 + ClassDefNV(ITS3Layer, 0); +}; +} // namespace o2::its3 #endif diff --git a/Detectors/Upgrades/ITS3/simulation/src/Detector.cxx b/Detectors/Upgrades/ITS3/simulation/src/Detector.cxx new file mode 100644 index 0000000000000..e69de29bb2d1d diff --git a/Detectors/Upgrades/ITS3/simulation/src/ITS3Layer.cxx b/Detectors/Upgrades/ITS3/simulation/src/ITS3Layer.cxx index 774b600ac3a09..4b107d9daa7db 100644 --- a/Detectors/Upgrades/ITS3/simulation/src/ITS3Layer.cxx +++ b/Detectors/Upgrades/ITS3/simulation/src/ITS3Layer.cxx @@ -12,378 +12,313 @@ /// \file ITS3Layer.h /// \brief Definition of the ITS3Layer class /// \author Fabrizio Grosa +/// \author felix.schlepper@cern.ch +#include "TGeoTube.h" +#include "TGeoVolume.h" + +#include "CommonConstants/MathConstants.h" +#include "ITSBase/Specs.h" #include "ITSBase/GeometryTGeo.h" #include "ITS3Base/SegmentationSuperAlpide.h" #include "ITS3Simulation/ITS3Layer.h" +#include "fairlogger/Logger.h" -#include // for LOG - -#include // for TGeoTube, TGeoTubeSeg -#include -#include // for TGeoVolume, TGeoVolumeAssembly - -using namespace o2::its3; - -/// \cond CLASSIMP -ClassImp(ITS3Layer); -/// \endcond +namespace o2::its3 +{ -ITS3Layer::ITS3Layer(int lay) - : TObject(), - mLayerNumber(lay) +void ITS3Layer::init() { - SegmentationSuperAlpide seg(lay); - mChipThickness = seg.mDetectorLayerThickness; - mSensorThickness = seg.mSensorLayerThickness; + // First we start by creating variables which we are reusing a couple of times. + mR = constants::radii[mNLayer]; + mRmin = mR - constants::thickness / 2.; + mRmax = mR + constants::thickness / 2.; } -ITS3Layer::~ITS3Layer() = default; - -void ITS3Layer::createLayer(TGeoVolume* motherVolume, double radiusBetweenLayer) +void ITS3Layer::createLayer(TGeoVolume* motherVolume, int layer) { - TGeoMedium* medSi = gGeoManager->GetMedium("IT3_SI$"); - TGeoMedium* medAir = gGeoManager->GetMedium("IT3_AIR$"); - - double rmin = mRadius; - double rmax = rmin + mChipThickness; - double rminSensor = rmax - mSensorThickness; - radiusBetweenLayer = radiusBetweenLayer - mChipThickness; - double phiGap = TMath::ASin(mGapPhi / 2.f / mRadius) * TMath::RadToDeg(); // degrees - double piDeg = TMath::Pi() * TMath::RadToDeg(); // degrees - - const int nElements = 7; - std::string names[nElements]; - - // do we need to keep the hierarchy? - names[0] = Form("%s%d", o2::its::GeometryTGeo::getITS3SensorPattern(), mLayerNumber); - names[1] = Form("%s%d", o2::its::GeometryTGeo::getITS3ChipPattern(), mLayerNumber); - names[2] = Form("%s%d", o2::its::GeometryTGeo::getITS3ModulePattern(), mLayerNumber); - names[3] = Form("%s%d", o2::its::GeometryTGeo::getITS3HalfStavePattern(), mLayerNumber); - names[4] = Form("%s%d", o2::its::GeometryTGeo::getITS3StavePattern(), mLayerNumber); - names[5] = Form("%s%d", o2::its::GeometryTGeo::getITS3HalfBarrelPattern(), mLayerNumber); - names[6] = Form("%s%d", o2::its::GeometryTGeo::getITS3LayerPattern(), mLayerNumber); - - TGeoTubeSeg* halfLayer[nElements - 1]; - TGeoVolume* volHalfLayer[nElements - 1]; - for (int iEl{0}; iEl < nElements - 1; ++iEl) { - TGeoMedium* med = (iEl <= 2) ? medSi : medAir; - if (iEl == 3) { - halfLayer[iEl] = new TGeoTubeSeg(rmin, rmax + radiusBetweenLayer, mZLen / 2, phiGap, piDeg - phiGap); - volHalfLayer[iEl] = new TGeoVolume(names[iEl].data(), halfLayer[iEl], med); - createCarbonFoamStructure(volHalfLayer[iEl], radiusBetweenLayer); - volHalfLayer[iEl]->SetVisibility(true); - volHalfLayer[iEl]->SetLineColor(kGray + 2); - } else if (iEl < 2) { - halfLayer[iEl] = new TGeoTubeSeg(rminSensor, rmax, mZLen / 2, phiGap, piDeg - phiGap); - volHalfLayer[iEl] = new TGeoVolume(names[iEl].data(), halfLayer[iEl], med); - volHalfLayer[iEl]->SetVisibility(true); - volHalfLayer[iEl]->SetLineColor(kRed + 1); - } else if (iEl == 2) { - halfLayer[iEl] = new TGeoTubeSeg(rmin, rmax, mZLen / 2, phiGap, piDeg - phiGap); - volHalfLayer[iEl] = new TGeoVolume(names[iEl].data(), halfLayer[iEl], med); - volHalfLayer[iEl]->SetVisibility(true); - volHalfLayer[iEl]->SetLineColor(kRed + 1); - } else { // all the others are simply half cylinders filling all the space - halfLayer[iEl] = new TGeoTubeSeg(rmin, rmax + radiusBetweenLayer, mZLen / 2, phiGap, piDeg - phiGap); - volHalfLayer[iEl] = new TGeoVolume(names[iEl].data(), halfLayer[iEl], med); - } - volHalfLayer[iEl]->SetUniqueID(mChipTypeID); - - if (iEl > 0) { - LOGP(debug, "Inserting {} inside {}", volHalfLayer[iEl - 1]->GetName(), volHalfLayer[iEl]->GetName()); - int id = 0; - if (iEl == 1) { - id = 1; - } - volHalfLayer[iEl]->AddNode(volHalfLayer[iEl - 1], id, nullptr); - } + // Create one layer of ITS3 and attach it to the motherVolume. + if (layer > 2) { + LOGP(fatal, "Cannot create more than 3 layers!"); + return; } + mNLayer = layer; + + init(); + createPixelArray(); + createTile(); + createRSU(); + createSegment(); + createChip(); + createCarbonForm(); + createLayerImpl(); + + // Add it to motherVolume + motherVolume->AddNode(mLayer, 0); +} - TGeoTranslation* translationTop = new TGeoTranslation(0., mGapY / 2, 0.); - TGeoTranslation* translationBottom = new TGeoTranslation(0., -mGapY / 2, 0.); - TGeoRotation* rotationBottom = new TGeoRotation("", piDeg, 0., 0.); - TGeoCombiTrans* rotoTranslationBottom = - new TGeoCombiTrans(*translationBottom, *rotationBottom); - - TGeoVolumeAssembly* volLayer = new TGeoVolumeAssembly(names[nElements - 1].data()); - volLayer->AddNode(volHalfLayer[nElements - 2], 0, translationTop); - volLayer->AddNode(volHalfLayer[nElements - 2], 1, rotoTranslationBottom); - - // Finally put everything in the mother volume - LOGP(debug, "Inserting {} inside {}", volLayer->GetName(), motherVolume->GetName()); - motherVolume->AddNode(volLayer, 1, nullptr); +void ITS3Layer::createPixelArray() +{ + using namespace constants::pixelarray; + // A pixel array is pure silicon and the sensitive part of our detector. + // It will be segmented into a 440x144 matrix by the + // SuperSegmentationAlpide. + // Pixel Array is just a longer version of the biasing but starts in phi at + // biasPhi2. + double pixelArrayPhi1 = + constants::tile::biasing::length / mRmin * constants::math::Rad2Deg; + double pixelArrayPhi2 = + length / mRmin * constants::math::Rad2Deg + pixelArrayPhi1; + auto pixelArray = new TGeoTubeSeg(mRmin, mRmax, width / 2., + pixelArrayPhi1, pixelArrayPhi2); + mPixelArray = new TGeoVolume( + Form("pixelarray_%d", mNLayer) /* TODO change to correct name */, + pixelArray, + gGeoManager->GetMedium(material::MaterialNames[material::Silicon])); + mPixelArray->SetLineColor(color); + mPixelArray->RegisterYourself(); } -void ITS3Layer::createLayerWithDeadZones(TGeoVolume* motherVolume, double radiusBetweenLayer) +void ITS3Layer::createTile() { - TGeoMedium* medSi = gGeoManager->GetMedium("IT3_SI$"); - TGeoMedium* medAir = gGeoManager->GetMedium("IT3_AIR$"); - - double rmin = mRadius; - double rmax = rmin + mChipThickness; - double rminSensor = rmax - mSensorThickness; - double rmed = (rmax + rmin) / 2; - // width of sensors of layers is calculated from r and chips' widths - double widthSensor = (TMath::Pi() * rmed - (mNumSubSensorsHalfLayer - 1) * mMiddleChipWidth - 2 * mFringeChipWidth) / mNumSubSensorsHalfLayer; - radiusBetweenLayer = radiusBetweenLayer - mChipThickness; - double phiGap = TMath::ASin(mGapPhi / 2.f / mRadius) * TMath::RadToDeg(); // degrees - double piDeg = TMath::Pi() * TMath::RadToDeg(); // degrees - - const int nElements = 7; - std::string names[nElements]; - int nObjPerElement[nElements] = {mNumSubSensorsHalfLayer, 1, 1, 1, 1, 1, 1}; // mNumSubSensorsHalfLayer chips and sensors per half layer - - names[0] = Form("%s%d", o2::its::GeometryTGeo::getITS3SensorPattern(), mLayerNumber); - names[1] = Form("%s%d", o2::its::GeometryTGeo::getITS3ChipPattern(), mLayerNumber); - names[2] = Form("%s%d", o2::its::GeometryTGeo::getITS3ModulePattern(), mLayerNumber); - names[3] = Form("%s%d", o2::its::GeometryTGeo::getITS3HalfStavePattern(), mLayerNumber); - names[4] = Form("%s%d", o2::its::GeometryTGeo::getITS3StavePattern(), mLayerNumber); - names[5] = Form("%s%d", o2::its::GeometryTGeo::getITS3HalfBarrelPattern(), mLayerNumber); - names[6] = Form("%s%d", o2::its::GeometryTGeo::getITS3LayerPattern(), mLayerNumber); - - std::array, nElements - 1> halfLayer{}; - TGeoVolume* volHalfLayer[nElements - 1]; - - for (int iEl{0}; iEl < nElements - 1; ++iEl) { - TGeoMedium* med = (iEl <= 2) ? medSi : medAir; - - for (int iObj{0}; iObj < nObjPerElement[iEl]; ++iObj) { - if (iEl == 0) { // subsensors (mNumSubSensorsHalfLayer sectors with dead zones) - if (iObj == 0) { - halfLayer[iEl].push_back(new TGeoTubeSeg(Form("subsens%dlayer%d", iObj, mLayerNumber), rminSensor, rmax, mZLen / 2, TMath::RadToDeg() * mFringeChipWidth / rmed + phiGap, TMath::RadToDeg() * (mFringeChipWidth + widthSensor) / rmed)); - } else if (iObj == mNumSubSensorsHalfLayer - 1) { - halfLayer[iEl].push_back(new TGeoTubeSeg(Form("subsens%dlayer%d", iObj, mLayerNumber), rminSensor, rmax, mZLen / 2, TMath::RadToDeg() * (mFringeChipWidth + iObj * widthSensor + iObj * mMiddleChipWidth) / rmed, piDeg - TMath::RadToDeg() * mFringeChipWidth / rmed - phiGap)); - } else { - halfLayer[iEl].push_back(new TGeoTubeSeg(Form("subsens%dlayer%d", iObj, mLayerNumber), rminSensor, rmax, mZLen / 2, TMath::RadToDeg() * (mFringeChipWidth + iObj * widthSensor + iObj * mMiddleChipWidth) / rmed, TMath::RadToDeg() * (mFringeChipWidth + (iObj + 1) * widthSensor + iObj * mMiddleChipWidth) / rmed)); - } - } else if (iEl == 1) { - halfLayer[iEl].push_back(new TGeoTubeSeg(rminSensor, rmax, mZLen / 2, phiGap, piDeg - phiGap)); - } else if (iEl == 2) { - halfLayer[iEl].push_back(new TGeoTubeSeg(rmin, rmax, mZLen / 2, phiGap, piDeg - phiGap)); - } else { // all the others are simply half cylinders filling all the space - halfLayer[iEl].push_back(new TGeoTubeSeg(rmin, rmax + radiusBetweenLayer, mZLen / 2, phiGap, piDeg - phiGap)); - } - } - - if (iEl == 0) { - std::string subSensNames = ""; - for (int iObj{0}; iObj < nObjPerElement[iEl] - 1; ++iObj) { - subSensNames += Form("subsens%dlayer%d+", iObj, mLayerNumber); - } - subSensNames += Form("subsens%dlayer%d", nObjPerElement[iEl] - 1, mLayerNumber); - TGeoCompositeShape* sensor = new TGeoCompositeShape(subSensNames.data()); - volHalfLayer[iEl] = new TGeoVolume(names[iEl].data(), sensor, med); - volHalfLayer[iEl]->SetUniqueID(mChipTypeID); - volHalfLayer[iEl]->SetVisibility(true); - volHalfLayer[iEl]->SetLineColor(kRed + 1); - } else { - volHalfLayer[iEl] = new TGeoVolume(names[iEl].data(), halfLayer[iEl][0], med); - volHalfLayer[iEl]->SetUniqueID(mChipTypeID); - if (iEl == 3) { - createCarbonFoamStructure(volHalfLayer[iEl], radiusBetweenLayer); - volHalfLayer[iEl]->SetVisibility(true); - volHalfLayer[iEl]->SetLineColor(kGray + 2); - } else if (iEl == 1) { - volHalfLayer[iEl]->SetVisibility(true); - volHalfLayer[iEl]->SetLineColor(kBlue + 2); - } - - int id = (iEl == 1) ? 1 : 0; - LOGP(debug, "Inserting {} id {} inside {}", volHalfLayer[iEl - 1]->GetName(), id, volHalfLayer[iEl]->GetName()); - volHalfLayer[iEl]->AddNode(volHalfLayer[iEl - 1], id, nullptr); - } + using namespace constants::tile; + // This functions creates a single Tile, which is the basic building block + // of the chip. It consists of a pixelArray (sensitive area), biasing, power + // switches and readout periphery (latter three are insensitive). We start + // building the tile with the left upper edge of the biasing as center of + // the tile’s z-coordinate axis. + mTile = new TGeoVolumeAssembly(Form("tile_%d", mNLayer)); + mTile->VisibleDaughters(); + + // Biasing + auto zMoveBiasing = new TGeoTranslation(0, 0, +biasing::width / 2.); + double biasPhi1 = 0; + double biasPhi2 = biasing::length / mRmin * constants::math::Rad2Deg; + auto biasing = + new TGeoTubeSeg(mRmin, mRmax, biasing::width / 2, biasPhi1, + biasPhi2); + auto biasingVol = new TGeoVolume( + Form("biasing_%d", mNLayer), biasing, + gGeoManager->GetMedium(material::MaterialNames[material::DeadZone])); + biasingVol->SetLineColor(biasing::color); + biasingVol->RegisterYourself(); + if (mVerbose) { + std::cout << "Biasing:" << std::endl; + biasingVol->InspectShape(); + biasingVol->InspectMaterial(); } + mTile->AddNode(biasingVol, 0, zMoveBiasing); + + // Pixel Array is just a longer version of the biasing but starts in phi at + // biasPhi2. + mTile->AddNode(mPixelArray, 0, zMoveBiasing); + + // The readout periphery is also on top of the pixel array but extrudes on +z a bit e.g. is wider. + auto zMoveReadout = new TGeoTranslation(0, 0, +readout::width / 2.); + double readoutPhi1 = + constants::pixelarray::length / mRmin * constants::math::Rad2Deg + biasPhi2; + double readoutPhi2 = + readout::length / mRmin * constants::math::Rad2Deg + readoutPhi1; + auto readout = new TGeoTubeSeg(mRmin, mRmax, readout::width / 2, + readoutPhi1, readoutPhi2); + auto readoutVol = new TGeoVolume( + Form("readout_%d", mNLayer), readout, + gGeoManager->GetMedium(material::MaterialNames[material::DeadZone])); + readoutVol->SetLineColor(readout::color); + readoutVol->RegisterYourself(); + if (mVerbose) { + std::cout << "Readout:" << std::endl; + readoutVol->InspectShape(); + readoutVol->InspectMaterial(); + } + mTile->AddNode(readoutVol, 0, zMoveReadout); + + // Power Switches are on the side right side of the pixel array and biasing. + auto zMovePowerSwitches = new TGeoTranslation(0, 0, +powerswitches::width / 2. + biasing::width); + double powerPhi1 = 0; + double powerPhi2 = powerswitches::length / mRmin * constants::math::Rad2Deg; + auto powerSwitches = new TGeoTubeSeg( + mRmin, mRmax, powerswitches::width / 2, powerPhi1, powerPhi2); + auto powerSwitchesVol = new TGeoVolume( + Form("powerswitches_%d", mNLayer), powerSwitches, + gGeoManager->GetMedium(material::MaterialNames[material::DeadZone])); + powerSwitchesVol->SetLineColor(powerswitches::color); + if (mVerbose) { + std::cout << "PowerSwitches:" << std::endl; + powerSwitchesVol->InspectShape(); + powerSwitchesVol->InspectMaterial(); + } + mTile->AddNode(powerSwitchesVol, 0, zMovePowerSwitches); - TGeoTranslation* translationTop = new TGeoTranslation(0., mGapY / 2, 0.); - TGeoTranslation* translationBottom = new TGeoTranslation(0., -mGapY / 2, 0.); - TGeoRotation* rotationBottom = new TGeoRotation("", piDeg, 0., 0.); - TGeoCombiTrans* rotoTranslationBottom = - new TGeoCombiTrans(*translationBottom, *rotationBottom); - - TGeoVolumeAssembly* volLayer = new TGeoVolumeAssembly(names[nElements - 1].data()); - volLayer->AddNode(volHalfLayer[nElements - 2], 0, translationTop); - volLayer->AddNode(volHalfLayer[nElements - 2], 1, rotoTranslationBottom); - - // Finally put everything in the mother volume - LOGP(debug, "Inserting {} inside {}", volLayer->GetName(), motherVolume->GetName()); - motherVolume->AddNode(volLayer, 1, nullptr); + if (mSubstrate) { + // Create the substrate layer at the back of the tile. + // TODO + } } -void ITS3Layer::create4thLayer(TGeoVolume* motherVolume) +void ITS3Layer::createRSU() { - TGeoMedium* medSi = gGeoManager->GetMedium("IT3_SI$"); - TGeoMedium* medAir = gGeoManager->GetMedium("IT3_AIR$"); - - double rmin = mRadius; - double rmax = rmin + mChipThickness; - double rminSensor = rmax - mSensorThickness; - double rmed = (rmax + rmin) / 2; - // width of sensors of layers is calculated from r and chips' widths - double widthSensor = (0.5 * TMath::Pi() * rmed - (mNumSubSensorsHalfLayer - 1) * mMiddleChipWidth - 2 * mFringeChipWidth) / mNumSubSensorsHalfLayer; - double phiGap = TMath::ASin(mGapPhi / 2.f / mRadius) * TMath::RadToDeg(); // degrees - double piDeg = TMath::Pi() * TMath::RadToDeg(); // degrees - - const int nElements = 7; - std::string names[nElements]; - int nObjPerElement[nElements] = {mNumSubSensorsHalfLayer, 1, 1, 1, 1, 1, 1}; // mNumSubSensorsHalfLayer chips and sensors per half layer - - names[0] = Form("%s%d", o2::its::GeometryTGeo::getITS3SensorPattern(), mLayerNumber); - names[1] = Form("%s%d", o2::its::GeometryTGeo::getITS3ChipPattern(), mLayerNumber); - names[2] = Form("%s%d", o2::its::GeometryTGeo::getITS3ModulePattern(), mLayerNumber); - names[3] = Form("%s%d", o2::its::GeometryTGeo::getITS3HalfStavePattern(), mLayerNumber); - names[4] = Form("%s%d", o2::its::GeometryTGeo::getITS3StavePattern(), mLayerNumber); - names[5] = Form("%s%d", o2::its::GeometryTGeo::getITS3HalfBarrelPattern(), mLayerNumber); - names[6] = Form("%s%d", o2::its::GeometryTGeo::getITS3LayerPattern(), mLayerNumber); - - std::array, nElements - 1> quarterLayer{}; - TGeoVolume* volQuarterLayer[nElements - 1]; - - for (int iEl{0}; iEl < nElements - 2; ++iEl) { // we need only 2 half barrels as usual, so we stop at the stave - TGeoMedium* med = (iEl <= 2) ? medSi : medAir; - - for (int iObj{0}; iObj < nObjPerElement[iEl]; ++iObj) { - if (iEl == 0) { // subsensors (mNumSubSensorsHalfLayer sectors with dead zones) - if (iObj == 0) { - quarterLayer[iEl].push_back(new TGeoTubeSeg(Form("subsens%dlayer%d", iObj, mLayerNumber), rminSensor, rmax, mZLen / 2, piDeg / 4. + TMath::RadToDeg() * mFringeChipWidth / rmed + phiGap, piDeg / 4. + TMath::RadToDeg() * (mFringeChipWidth + widthSensor) / rmed)); - } else if (iObj == mNumSubSensorsHalfLayer - 1) { - quarterLayer[iEl].push_back(new TGeoTubeSeg(Form("subsens%dlayer%d", iObj, mLayerNumber), rminSensor, rmax, mZLen / 2, piDeg / 4. + TMath::RadToDeg() * (mFringeChipWidth + iObj * widthSensor + iObj * mMiddleChipWidth) / rmed, 3. / 4. * piDeg - TMath::RadToDeg() * mFringeChipWidth / rmed - phiGap)); - } else { - quarterLayer[iEl].push_back(new TGeoTubeSeg(Form("subsens%dlayer%d", iObj, mLayerNumber), rminSensor, rmax, mZLen / 2, piDeg / 4. + TMath::RadToDeg() * (mFringeChipWidth + iObj * widthSensor + iObj * mMiddleChipWidth) / rmed, piDeg / 4. + TMath::RadToDeg() * (mFringeChipWidth + (iObj + 1) * widthSensor + iObj * mMiddleChipWidth) / rmed)); - } - } else if (iEl == 1) { - quarterLayer[iEl].push_back(new TGeoTubeSeg(rminSensor, rmax, mZLen / 2, phiGap, piDeg / 4. + piDeg * 3 / 4. - phiGap)); - } else if (iEl == 2) { - quarterLayer[iEl].push_back(new TGeoTubeSeg(rmin, rmax, mZLen / 2, piDeg / 4. + phiGap, piDeg * 3 / 4. - phiGap)); - } else { // all the others are simply quarter cylinders filling all the space - quarterLayer[iEl].push_back(new TGeoTubeSeg(rmin, rmax, mZLen / 2, piDeg / 4. + phiGap, piDeg * 3 / 4. - phiGap)); - } - } - - if (iEl == 0) { - std::string subSensNames = ""; - for (int iObj{0}; iObj < nObjPerElement[iEl] - 1; ++iObj) { - subSensNames += Form("subsens%dlayer%d+", iObj, mLayerNumber); - } - subSensNames += Form("subsens%dlayer%d", nObjPerElement[iEl] - 1, mLayerNumber); - TGeoCompositeShape* sensor = new TGeoCompositeShape(subSensNames.data()); - volQuarterLayer[iEl] = new TGeoVolume(names[iEl].data(), sensor, med); - volQuarterLayer[iEl]->SetUniqueID(mChipTypeID); - volQuarterLayer[iEl]->SetVisibility(true); - volQuarterLayer[iEl]->SetLineColor(kRed + 1); - } else { - volQuarterLayer[iEl] = new TGeoVolume(names[iEl].data(), quarterLayer[iEl][0], med); - volQuarterLayer[iEl]->SetUniqueID(mChipTypeID); - if (iEl == 4) { - // createCarbonFoamStructure(volQuarterLayer[iEl], 0., true); - volQuarterLayer[iEl]->SetVisibility(true); - volQuarterLayer[iEl]->SetLineColor(kGray + 2); - } else if (iEl == 1) { - volQuarterLayer[iEl]->SetVisibility(true); - volQuarterLayer[iEl]->SetLineColor(kBlue + 2); - } - - int id = (iEl == 1) ? 1 : 0; - LOGP(debug, "Inserting {} id {} inside {}", volQuarterLayer[iEl - 1]->GetName(), id, volQuarterLayer[iEl]->GetName()); - volQuarterLayer[iEl]->AddNode(volQuarterLayer[iEl - 1], id, nullptr); - } + using namespace constants::rsu; + // A Repeated Sensor Unit (RSU) is 12 Tiles + 4 Databackbones stichted together. + mRSU = new TGeoVolumeAssembly(Form("rsu_%d", mNLayer)); + mRSU->VisibleDaughters(); + int nCopyRSU{0}, nCopyDB{0}; + + // Create the DatabackBone + // The Databackbone spans the whole phi of the tile. + double dataBackbonePhi1 = 0; + double dataBackbonePhi2 = databackbone::length / mRmin * + constants::math::Rad2Deg; + auto dataBackbone = new TGeoTubeSeg(mRmin, mRmax, databackbone::width / 2., + dataBackbonePhi1, + dataBackbonePhi2); + auto dataBackboneVol = new TGeoVolume( + Form("databackbone_%d", mNLayer), dataBackbone, + gGeoManager->GetMedium(material::MaterialNames[material::DeadZone])); + dataBackboneVol->SetLineColor(databackbone::color); + dataBackboneVol->RegisterYourself(); + if (mVerbose) { + std::cout << "DataBackbone:" << std::endl; + dataBackboneVol->InspectShape(); + dataBackboneVol->InspectMaterial(); } - TGeoTranslation* translationTopRight = new TGeoTranslation(mGapXDirection / 2, 0., 0.); - TGeoRotation* rotationTopRight = new TGeoRotation("", -piDeg / 4, 0., 0.); - TGeoCombiTrans* rotoTranslationTopRight = new TGeoCombiTrans(*translationTopRight, *rotationTopRight); - - TGeoTranslation* translationTopLeft = new TGeoTranslation(-mGapXDirection / 2, 0., 0.); - TGeoRotation* rotationTopLeft = new TGeoRotation("", piDeg / 4, 0., 0.); - TGeoCombiTrans* rotoTranslationTopLeft = new TGeoCombiTrans(*translationTopLeft, *rotationTopLeft); - - TGeoVolumeAssembly* halfLayer = new TGeoVolumeAssembly(names[nElements - 2].data()); - halfLayer->AddNode(volQuarterLayer[nElements - 3], 0, rotoTranslationTopRight); - halfLayer->AddNode(volQuarterLayer[nElements - 3], 1, rotoTranslationTopLeft); - - TGeoTranslation* translationTop = new TGeoTranslation(0., mGapY / 2, 0.); - TGeoTranslation* translationBottom = new TGeoTranslation(0., -mGapY / 2, 0.); - TGeoRotation* rotationBottom = new TGeoRotation("", piDeg, 0., 0.); - TGeoCombiTrans* rotoTranslationBottom = new TGeoCombiTrans(*translationBottom, *rotationBottom); - - TGeoVolumeAssembly* volLayer = new TGeoVolumeAssembly(names[nElements - 1].data()); - volLayer->AddNode(halfLayer, 0, translationTop); - volLayer->AddNode(halfLayer, 1, rotoTranslationBottom); - - // Finally put everything in the mother volume - LOGP(debug, "Inserting {} inside {}", volLayer->GetName(), motherVolume->GetName()); - motherVolume->AddNode(volLayer, 1, nullptr); + // Lower Left + auto zMoveLL1 = new TGeoTranslation(0, 0, constants::tile::width); + auto zMoveLL2 = new TGeoTranslation(0, 0, constants::tile::width * 2.); + auto zMoveLLDB = new TGeoTranslation(0, 0, -databackbone::width / 2.); + // Lets attach the tiles to the QS. + mRSU->AddNode(mTile, nCopyRSU++, nullptr); + mRSU->AddNode(mTile, nCopyRSU++, zMoveLL1); + mRSU->AddNode(mTile, nCopyRSU++, zMoveLL2); + mRSU->AddNode(dataBackboneVol, nCopyDB++, zMoveLLDB); + + // Lower Right + auto zMoveLR0 = new TGeoTranslation(0, 0, +width / 2.); + auto zMoveLR1 = new TGeoTranslation(0, 0, constants::tile::width + width / 2.); + auto zMoveLR2 = new TGeoTranslation(0, 0, constants::tile::width * 2. + width / 2.); + auto zMoveLRDB = new TGeoTranslation(0, 0, -databackbone::width / 2. + width / 2.); + // Lets attach the tiles to the QS. + mRSU->AddNode(mTile, nCopyRSU++, zMoveLR0); + mRSU->AddNode(mTile, nCopyRSU++, zMoveLR1); + mRSU->AddNode(mTile, nCopyRSU++, zMoveLR2); + mRSU->AddNode(dataBackboneVol, nCopyDB++, zMoveLRDB); + + // Rotation for top half + double phi = length / mRmin * constants::math::Rad2Deg; + auto rot = new TGeoRotation("", 0, 0, phi / 2.); + + // Upper Left + auto zMoveUL1 = new TGeoCombiTrans(0, 0, constants::tile::width, rot); + auto zMoveUL2 = new TGeoCombiTrans(0, 0, constants::tile::width * 2., rot); + auto zMoveULDB = new TGeoCombiTrans(0, 0, -databackbone::width / 2., rot); + // Lets attach the tiles to the QS. + mRSU->AddNode(mTile, nCopyRSU++, rot); + mRSU->AddNode(mTile, nCopyRSU++, zMoveUL1); + mRSU->AddNode(mTile, nCopyRSU++, zMoveUL2); + mRSU->AddNode(dataBackboneVol, nCopyDB++, zMoveULDB); + + // Upper Right + auto zMoveUR0 = new TGeoCombiTrans(0, 0, +width / 2., rot); + auto zMoveUR1 = new TGeoCombiTrans(0, 0, constants::tile::width + width / 2., rot); + auto zMoveUR2 = new TGeoCombiTrans(0, 0, constants::tile::width * 2. + width / 2., rot); + auto zMoveURDB = new TGeoCombiTrans(0, 0, -databackbone::width / 2. + width / 2., rot); + // Lets attach the tiles to the QS. + mRSU->AddNode(mTile, nCopyRSU++, zMoveUR0); + mRSU->AddNode(mTile, nCopyRSU++, zMoveUR1); + mRSU->AddNode(mTile, nCopyRSU++, zMoveUR2); + mRSU->AddNode(dataBackboneVol, nCopyDB++, zMoveURDB); } -void ITS3Layer::createCarbonFoamStructure(TGeoVolume* motherVolume, double deltaR, bool fourthLayer) +void ITS3Layer::createSegment() { - TGeoMedium* medCarbonFoam = (mBuildLevel < 1) ? gGeoManager->GetMedium("IT3_ERGDUOCEL$") : gGeoManager->GetMedium("IT3_AIR$"); // if build level >= 1 we do not put carbon foam but air - TGeoMedium* medGlue = (mBuildLevel < 2) ? gGeoManager->GetMedium("IT3_IMPREG_FLEECE$") : gGeoManager->GetMedium("IT3_AIR$"); // if build level >= 2 we do not put glue but air - TGeoMedium* medSi = (mBuildLevel <= 2) ? gGeoManager->GetMedium("IT3_SI$") : gGeoManager->GetMedium("IT3_AIR$"); // if build level > 2 we do not put silicon but air - - double rmaxWoAddMat = mRadius + mChipThickness; - double rmax = rmaxWoAddMat + mAddMaterial; - double radiusBetweenLayer = deltaR - mChipThickness; - double rmedFoam = rmax + radiusBetweenLayer / 2; - double phiGap = TMath::ASin(mGapPhi / 2.f / mRadius) * TMath::RadToDeg(); // degrees - double piDeg = TMath::Pi() * TMath::RadToDeg(); // degrees - double phiMin = (fourthLayer) ? piDeg / 4 : 0.; // degrees - double phiMax = (fourthLayer) ? piDeg * 3 / 4 : piDeg; // degrees - - TGeoTranslation* transSemicircle[2]; - transSemicircle[0] = new TGeoTranslation("transSemicircleFoam0", 0, 0, (mZLen - mLengthSemiCircleFoam) / 2); - transSemicircle[0]->RegisterYourself(); - transSemicircle[1] = new TGeoTranslation("transSemicircleFoam1", 0, 0, -(mZLen - mLengthSemiCircleFoam) / 2); - transSemicircle[1]->RegisterYourself(); - - TGeoTubeSeg* subGluedFoamBottom[4]; - subGluedFoamBottom[0] = new TGeoTubeSeg(Form("subgluedfoambottom0layer%d", mLayerNumber), rmax, rmax + mThickGluedFoam, mLengthSemiCircleFoam / 2, phiMin + phiGap, phiMax - phiGap); - subGluedFoamBottom[1] = new TGeoTubeSeg(Form("subgluedfoambottom1layer%d", mLayerNumber), rmax, rmax + mThickGluedFoam, mLengthSemiCircleFoam / 2, phiMin + phiGap, phiMax - phiGap); - subGluedFoamBottom[2] = new TGeoTubeSeg(Form("subgluedfoambottom2layer%d", mLayerNumber), rmax, rmax + mThickGluedFoam, (mZLen - mLengthSemiCircleFoam) / 2, phiMin + phiGap, phiMin + TMath::RadToDeg() * mHeightStripFoam / rmedFoam + phiGap); - subGluedFoamBottom[3] = new TGeoTubeSeg(Form("subgluedfoambottom3layer%d", mLayerNumber), rmax, rmax + mThickGluedFoam, (mZLen - mLengthSemiCircleFoam) / 2, phiMax - TMath::RadToDeg() * (mHeightStripFoam / rmedFoam) - phiGap, phiMax - phiGap); - TGeoTubeSeg* subGluedFoamTop[4]; - subGluedFoamTop[0] = new TGeoTubeSeg(Form("subgluedfoamtop0layer%d", mLayerNumber), rmax + radiusBetweenLayer - mThickGluedFoam, rmax + radiusBetweenLayer, mLengthSemiCircleFoam / 2, phiMin + phiGap, phiMax - phiGap); - subGluedFoamTop[1] = new TGeoTubeSeg(Form("subgluedfoamtop1layer%d", mLayerNumber), rmax + radiusBetweenLayer - mThickGluedFoam, rmax + radiusBetweenLayer, mLengthSemiCircleFoam / 2, phiMin + phiGap, phiMax - phiGap); - subGluedFoamTop[2] = new TGeoTubeSeg(Form("subgluedfoamtop2layer%d", mLayerNumber), rmax + radiusBetweenLayer - mThickGluedFoam, rmax + radiusBetweenLayer, (mZLen - mLengthSemiCircleFoam) / 2, phiMin + phiGap, phiMin + TMath::RadToDeg() * mHeightStripFoam / rmedFoam + phiGap); - subGluedFoamTop[3] = new TGeoTubeSeg(Form("subgluedfoamtop3layer%d", mLayerNumber), rmax + radiusBetweenLayer - mThickGluedFoam, rmax + radiusBetweenLayer, (mZLen - mLengthSemiCircleFoam) / 2, phiMax - TMath::RadToDeg() * (mHeightStripFoam / rmedFoam) - phiGap, phiMax - phiGap); - - std::string subGluedFoamsNames = ""; - for (int iObj{0}; iObj < 2; ++iObj) { - subGluedFoamsNames += Form("(subgluedfoambottom%dlayer%d:transSemicircleFoam%d)+", iObj, mLayerNumber, iObj); + using namespace constants::segment; + // A segment is 12 RSUs + left and right end cap. We place the first rsu + // as z-coordinate center and attach to this. Hence, we will displace the + // left end-cap to the left and the right to right. + mSegment = new TGeoVolumeAssembly(Form("segment_%d", mNLayer)); + mSegment->VisibleDaughters(); + + for (int i{0}; i < nRSUs; ++i) { + auto zMove = new TGeoTranslation(0, 0, +i * constants::rsu::width + constants::rsu::databackbone::width); + mSegment->AddNode(mRSU, i, zMove); } - subGluedFoamsNames += Form("subgluedfoambottom2layer%d+", mLayerNumber); - subGluedFoamsNames += Form("subgluedfoambottom3layer%d+", mLayerNumber); - for (int iObj{0}; iObj < 2; ++iObj) { - subGluedFoamsNames += Form("(subgluedfoamtop%dlayer%d:transSemicircleFoam%d)+", iObj, mLayerNumber, iObj); + // LEC + double lecPhi1 = 0; + double lecPhi2 = lec::length / mRmin * constants::math::Rad2Deg; + auto zMoveLEC = new TGeoTranslation(0, 0, -lec::width / 2.); + auto lec = + new TGeoTubeSeg(mRmin, mRmax, lec::width / 2., lecPhi1, lecPhi2); + auto lecVol = new TGeoVolume( + Form("lec_%d", mNLayer), lec, + gGeoManager->GetMedium(material::MaterialNames[material::DeadZone])); + lecVol->SetLineColor(lec::color); + lecVol->RegisterYourself(); + if (mVerbose) { + std::cout << "LEC:" << std::endl; + lecVol->InspectShape(); + lecVol->InspectMaterial(); } - subGluedFoamsNames += Form("subgluedfoamtop2layer%d+", mLayerNumber); - subGluedFoamsNames += Form("subgluedfoamtop3layer%d", mLayerNumber); - - TGeoCompositeShape* gluedfoam = new TGeoCompositeShape(subGluedFoamsNames.data()); - TGeoVolume* volGlue = new TGeoVolume(Form("Glue%d", mLayerNumber), gluedfoam, medGlue); - motherVolume->AddNode(volGlue, 1, nullptr); - - TGeoTubeSeg* subFoam[4]; - subFoam[0] = new TGeoTubeSeg(Form("subfoam0layer%d", mLayerNumber), rmax + mThickGluedFoam, rmax + radiusBetweenLayer - mThickGluedFoam, mLengthSemiCircleFoam / 2, phiMin + phiGap, phiMax - phiGap); - subFoam[1] = new TGeoTubeSeg(Form("subfoam1layer%d", mLayerNumber), rmax + mThickGluedFoam, rmax + radiusBetweenLayer - mThickGluedFoam, mLengthSemiCircleFoam / 2, phiMin + phiGap, phiMax - phiGap); - subFoam[2] = new TGeoTubeSeg(Form("subfoam2layer%d", mLayerNumber), rmax + mThickGluedFoam, rmax + radiusBetweenLayer - mThickGluedFoam, (mZLen - mLengthSemiCircleFoam) / 2, phiMin + phiGap, phiMin + TMath::RadToDeg() * mHeightStripFoam / rmedFoam + phiGap); - subFoam[3] = new TGeoTubeSeg(Form("subfoam3layer%d", mLayerNumber), rmax + mThickGluedFoam, rmax + radiusBetweenLayer - mThickGluedFoam, (mZLen - mLengthSemiCircleFoam) / 2, phiMax - TMath::RadToDeg() * (mHeightStripFoam / rmedFoam) - phiGap, phiMax - phiGap); - - std::string subFoamNames = ""; - for (int iObj{0}; iObj < 2; ++iObj) { - subFoamNames += Form("(subfoam%dlayer%d:transSemicircleFoam%d)+", iObj, mLayerNumber, iObj); + mSegment->AddNode(lecVol, 0, zMoveLEC); + + // REC; reuses lecPhi1,2 + auto zMoveREC = new TGeoTranslation(0, 0, nRSUs * constants::rsu::width + rec::width / 2.); + auto rec = + new TGeoTubeSeg(mRmin, mRmax, rec::width / 2., lecPhi1, lecPhi2); + auto recVol = new TGeoVolume( + Form("rec_%d", mNLayer), rec, + gGeoManager->GetMedium(material::MaterialNames[material::DeadZone])); + recVol->SetLineColor(rec::color); + recVol->RegisterYourself(); + if (mVerbose) { + std::cout << "REC:" << std::endl; + recVol->InspectShape(); + recVol->InspectMaterial(); } - subFoamNames += Form("subfoam2layer%d+", mLayerNumber); - subFoamNames += Form("subfoam3layer%d", mLayerNumber); + mSegment->AddNode(recVol, 0, zMoveREC); +} + +void ITS3Layer::createChip() +{ - TGeoCompositeShape* foam = new TGeoCompositeShape(subFoamNames.data()); - TGeoVolume* volFoam = new TGeoVolume(Form("CarbonFoam%d", mLayerNumber), foam, medCarbonFoam); - motherVolume->AddNode(volFoam, 1, nullptr); + // A HalfLayer is composed out of multiple segment stitched together along + // rphi. + mChip = new TGeoVolumeAssembly(Form("chip_%d", mNLayer)); + mChip->VisibleDaughters(); - if (mAddMaterial > 0.) { - TGeoTubeSeg* addMat = new TGeoTubeSeg(Form("additionalMaterialLayer%d", mLayerNumber), rmaxWoAddMat, rmax, mLengthSemiCircleFoam / 2, phiMin + phiGap, phiMax - phiGap); - TGeoVolume* volAddMat = new TGeoVolume(Form("AdditionalMaterial%d", mLayerNumber), addMat, medSi); - motherVolume->AddNode(volAddMat, 1, nullptr); + for (int i{0}; i < constants::nSegments[mNLayer]; ++i) { + double phiOffset = constants::segment::length / mRmin * constants::math::Rad2Deg; + auto rot = new TGeoRotation("", 0, 0, phiOffset * i); + mChip->AddNode(mSegment, i, rot); } } + +void ITS3Layer::createCarbonForm() +{ + mCarbonForm = new TGeoVolumeAssembly(Form("carbonform_%d", mNLayer)); + mCarbonForm->VisibleDaughters(); + mCarbonForm->AddNode(mChip, 0); + // TODO +} + +void ITS3Layer::createLayerImpl() +{ + // At long last a single layer... A layer is two HalfLayers (duuhhh) but + // we have to take care of the equatorial gap. So both half layers will be + // offset slightly by rotating in phi the upper HalfLayer and negative phi + // the other one. + mLayer = new TGeoVolumeAssembly(Form("layer_%d", mNLayer)); + mLayer->VisibleDaughters(); + + // The offset is the right angle triangle of the middle radius with the + // transverse axis. + double phiOffset = std::asin(constants::equatorialGap / mR) * constants::math::Rad2Deg; + // double phiOffset = constants::equatorialGap / mRmin / 2.; + auto rotTop = new TGeoRotation("", 0, 0, +phiOffset); + auto rotBot = new TGeoRotation("", 0, 0, phiOffset + 180); + + mLayer->AddNode(mCarbonForm, 0, rotTop); + mLayer->AddNode(mCarbonForm, 1, rotBot); +} +} // namespace o2::its3 diff --git a/Detectors/Vertexing/src/SVertexer.cxx b/Detectors/Vertexing/src/SVertexer.cxx index 3458e2c19f507..5454c93228330 100644 --- a/Detectors/Vertexing/src/SVertexer.cxx +++ b/Detectors/Vertexing/src/SVertexer.cxx @@ -626,7 +626,7 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP, if (mV0Hyps[Lambda].checkTight(p2Pos, p2Neg, p2V0, ptV0) && (!mSVParams->mRequireTPCforCascBaryons || seedP.hasTPC) && seedP.compatibleProton) { goodLamForCascade = true; } - if (mV0Hyps[AntiLambda].checkTight(p2Pos, p2Neg, p2V0, ptV0) && (!mSVParams->mRequireTPCforCascBaryons || seedN.hasTPC && seedN.compatibleProton)) { + if (mV0Hyps[AntiLambda].checkTight(p2Pos, p2Neg, p2V0, ptV0) && (!mSVParams->mRequireTPCforCascBaryons || seedN.hasTPC) && seedN.compatibleProton) { goodALamForCascade = true; } } diff --git a/Detectors/ZDC/calib/src/TDCCalibSpec.cxx b/Detectors/ZDC/calib/src/TDCCalibSpec.cxx index aa6cab936805d..248b18ccfdfc0 100644 --- a/Detectors/ZDC/calib/src/TDCCalibSpec.cxx +++ b/Detectors/ZDC/calib/src/TDCCalibSpec.cxx @@ -199,8 +199,8 @@ framework::DataProcessorSpec getTDCCalibSpec() std::vector inputs; inputs.emplace_back("tdccalibconfig", "ZDC", "TDCCALIBCONFIG", 0, Lifetime::Condition, o2::framework::ccdbParamSpec(fmt::format("{}", o2::zdc::CCDBPathTDCCalibConfig.data()))); inputs.emplace_back("tdccalib", "ZDC", "TDCCALIB", 0, Lifetime::Condition, o2::framework::ccdbParamSpec(fmt::format("{}", o2::zdc::CCDBPathTDCCalib.data()))); - inputs.emplace_back("tdccalibdata", "ZDC", "TDCCALIBDATA", 0, Lifetime::Timeframe); - inputs.emplace_back("tdc_1dh", ConcreteDataTypeMatcher{"ZDC", "TDC_1DH"}, Lifetime::Timeframe); + inputs.emplace_back("tdccalibdata", "ZDC", "TDCCALIBDATA", 0, Lifetime::Sporadic); + inputs.emplace_back("tdc_1dh", ConcreteDataTypeMatcher{"ZDC", "TDC_1DH"}, Lifetime::Sporadic); std::vector outputs; outputs.emplace_back(ConcreteDataTypeMatcher{o2::calibration::Utils::gDataOriginCDBPayload, "ZDC_TDCcalib"}, Lifetime::Sporadic); diff --git a/Framework/Core/include/Framework/ASoA.h b/Framework/Core/include/Framework/ASoA.h index a0a45a833efeb..898273b63d043 100644 --- a/Framework/Core/include/Framework/ASoA.h +++ b/Framework/Core/include/Framework/ASoA.h @@ -614,14 +614,14 @@ struct DefaultIndexPolicy : IndexPolicyBase { this->setCursor(mMaxRow); } - bool operator!=(DefaultIndexPolicy const& other) const + friend bool operator!=(DefaultIndexPolicy const& lh, DefaultIndexPolicy const& rh) { - return O2_BUILTIN_LIKELY(this->mRowIndex != other.mRowIndex); + return O2_BUILTIN_LIKELY(lh.mRowIndex != rh.mRowIndex); } - bool operator==(DefaultIndexPolicy const& other) const + friend bool operator==(DefaultIndexPolicy const& lh, DefaultIndexPolicy const& rh) { - return O2_BUILTIN_UNLIKELY(this->mRowIndex == other.mRowIndex); + return O2_BUILTIN_UNLIKELY(lh.mRowIndex == rh.mRowIndex); } bool operator!=(RowViewSentinel const& sentinel) const diff --git a/Framework/Core/include/Framework/ASoAHelpers.h b/Framework/Core/include/Framework/ASoAHelpers.h index 1573ce3abf235..932f4fd029298 100644 --- a/Framework/Core/include/Framework/ASoAHelpers.h +++ b/Framework/Core/include/Framework/ASoAHelpers.h @@ -1253,13 +1253,13 @@ struct CombinationsGenerator { { return this->mCurrent; } - bool operator==(const CombinationsIterator& rh) + friend bool operator==(const CombinationsIterator& lh, const CombinationsIterator& rh) { - return (this->mIsEnd && rh.mIsEnd) || (this->mCurrent == rh.mCurrent); + return (lh.mIsEnd && rh.mIsEnd) || (lh.mCurrent == rh.mCurrent); } - bool operator!=(const CombinationsIterator& rh) + friend bool operator!=(const CombinationsIterator& lh, const CombinationsIterator& rh) { - return !(*this == rh); + return !(lh == rh); } }; diff --git a/Framework/Core/include/Framework/DataOutputDirector.h b/Framework/Core/include/Framework/DataOutputDirector.h index 86b23782a89a0..8661dc3eb87e5 100644 --- a/Framework/Core/include/Framework/DataOutputDirector.h +++ b/Framework/Core/include/Framework/DataOutputDirector.h @@ -13,9 +13,7 @@ #include "TFile.h" -#include "Framework/DataDescriptorQueryBuilder.h" #include "Framework/DataDescriptorMatcher.h" -#include "Framework/DataSpecUtils.h" #include "Framework/InputSpec.h" #include "rapidjson/fwd.h" diff --git a/Framework/Core/include/Framework/DataSpecUtils.h b/Framework/Core/include/Framework/DataSpecUtils.h index 3dcdf6d790cd9..65f8585302aa7 100644 --- a/Framework/Core/include/Framework/DataSpecUtils.h +++ b/Framework/Core/include/Framework/DataSpecUtils.h @@ -80,9 +80,24 @@ struct DataSpecUtils { /// @return true if the InputSpec will match at least the provided @a origin. static bool partialMatch(InputSpec const& spec, o2::header::DataOrigin const& origin); + /// @return true if the InputSpec will match at least one of the provided @a origins + template + static bool partialMatch(InputSpec const& spec, std::array const& origins) + { + return std::find_if(origins.begin(), origins.end(), [&](auto const& o) { return DataSpecUtils::asConcreteOrigin(spec) == o; }) != origins.end(); + } + /// @return true if the OutputSpec will match at least the provided @a origin. static bool partialMatch(OutputSpec const& spec, o2::header::DataOrigin const& origin); + /// @return true if the OutputSpec will match at least one of the provided @a origins + template + static bool partialMatch(OutputSpec const& spec, std::array const& origins) + { + auto dataType = DataSpecUtils::asConcreteDataTypeMatcher(spec); + return std::find_if(origins.begin(), origins.end(), [&](auto const& o) { return dataType.origin == o; }) != origins.end(); + } + /// @return true if the OutputSpec will match at least the provided @a description. static bool partialMatch(InputSpec const& spec, o2::header::DataDescription const& description); diff --git a/Framework/Core/include/Framework/GroupedCombinations.h b/Framework/Core/include/Framework/GroupedCombinations.h index 4d2d7678493dc..5904916613d2e 100644 --- a/Framework/Core/include/Framework/GroupedCombinations.h +++ b/Framework/Core/include/Framework/GroupedCombinations.h @@ -137,11 +137,11 @@ struct GroupedCombinationsGenerator { { return *mCurrentGrouped; } - bool operator==(const GroupedIterator& rh) + bool operator==(const GroupedIterator& rh) const { return (this->mIsEnd && rh.mIsEnd) || (this->mCurrent == rh.mCurrent); } - bool operator!=(const GroupedIterator& rh) + bool operator!=(const GroupedIterator& rh) const { return !(*this == rh); } diff --git a/Framework/Core/include/Framework/VariantJSONHelpers.h b/Framework/Core/include/Framework/VariantJSONHelpers.h index ce998275d26b3..ebf0b308581a1 100644 --- a/Framework/Core/include/Framework/VariantJSONHelpers.h +++ b/Framework/Core/include/Framework/VariantJSONHelpers.h @@ -201,7 +201,7 @@ struct VariantReader : public rapidjson::BaseReaderHandler, Va return false; } if (states.top() == State::IN_DATA) { - //no previous keys + // no previous keys states.push(State::IN_KEY); currentKey = str; return true; @@ -281,14 +281,14 @@ struct VariantReader : public rapidjson::BaseReaderHandler, Va return false; } if (states.top() == State::IN_ARRAY) { - //finish up array + // finish up array states.pop(); if constexpr (isArray2D() || isLabeledArray()) { rows = elementCount; } return true; } else if (states.top() == State::IN_ROW) { - //finish up row + // finish up row states.pop(); if constexpr (isArray2D() || isLabeledArray()) { cols = elementCount; @@ -350,6 +350,8 @@ void writeVariant(std::ostream& o, Variant const& v) w.Int(array2d(i, j)); } else if constexpr (std::is_same_v || std::is_same_v) { w.Double(array2d(i, j)); + } else if constexpr (std::is_same_v) { + w.String(array2d(i, j).c_str()); } } w.EndArray(); @@ -447,6 +449,6 @@ struct VariantJSONHelpers { } } }; -} +} // namespace o2::framework #endif // FRAMEWORK_VARIANTJSONHELPERS_H diff --git a/Framework/Core/src/CommonDataProcessors.cxx b/Framework/Core/src/CommonDataProcessors.cxx index 74e020595b345..e9cdc97e316ca 100644 --- a/Framework/Core/src/CommonDataProcessors.cxx +++ b/Framework/Core/src/CommonDataProcessors.cxx @@ -259,9 +259,11 @@ DataProcessorSpec CommonDataProcessors::getOutputObjHistSink(std::vector +#include #include "rapidjson/document.h" #include "rapidjson/prettywriter.h" @@ -532,7 +536,7 @@ bool DataOutputDirector::checkFileSizes() // loop over all files // if one file is large, then all files need to be closed - for (int i = 0; i < mfilenameBases.size(); i++) { + for (auto i = 0U; i < mfilenameBases.size(); i++) { if (!mfilePtrs[i]) { continue; } @@ -557,7 +561,7 @@ bool DataOutputDirector::checkFileSizes() void DataOutputDirector::closeDataFiles() { - for (int i = 0; i < mfilePtrs.size(); i++) { + for (auto i = 0U; i < mfilePtrs.size(); i++) { auto filePtr = mfilePtrs[i]; if (filePtr) { if (filePtr->IsOpen() && mParentMaps[i]->GetEntries() > 0) { @@ -633,6 +637,5 @@ void DataOutputDirector::setMaximumFileSize(float maxfs) { mmaxfilesize = maxfs; } - } // namespace framework } // namespace o2 diff --git a/Framework/Core/src/DataProcessingDevice.cxx b/Framework/Core/src/DataProcessingDevice.cxx index 4d0cac7087611..b03904c5b3776 100644 --- a/Framework/Core/src/DataProcessingDevice.cxx +++ b/Framework/Core/src/DataProcessingDevice.cxx @@ -2234,7 +2234,7 @@ bool DataProcessingDevice::tryDispatchComputation(ServiceRegistryRef ref, std::v buffer[ai] = record.isValid(ai) ? '2' : '0'; } buffer[record.size()] = 0; - states.updateState({.id = short((int)ProcessingStateId::DATA_RELAYER_BASE + action.slot.index), (int)(record.size() + buffer - relayerSlotState), relayerSlotState}); + states.updateState({.id = short((int)ProcessingStateId::DATA_RELAYER_BASE + action.slot.index), .size = (int)(record.size() + buffer - relayerSlotState), .data = relayerSlotState}); }; // This is the main dispatching loop diff --git a/Framework/Core/src/DataSpecUtils.cxx b/Framework/Core/src/DataSpecUtils.cxx index 7cc9b8d03b673..b0a20df065878 100644 --- a/Framework/Core/src/DataSpecUtils.cxx +++ b/Framework/Core/src/DataSpecUtils.cxx @@ -12,10 +12,10 @@ #include "Framework/DataDescriptorMatcher.h" #include "Framework/DataMatcherWalker.h" #include "Framework/VariantHelpers.h" -#include "Framework/Logger.h" #include "Framework/RuntimeError.h" #include "Headers/DataHeaderHelpers.h" +#include #include #include #include @@ -516,7 +516,7 @@ DataDescriptorMatcher DataSpecUtils::dataDescriptorMatcherFrom(ConcreteDataMatch SubSpecificationTypeValueMatcher{concrete.subSpec}, std::make_unique(DataDescriptorMatcher::Op::Just, StartTimeValueMatcher{ContextRef{0}})))}; - return std::move(matchEverything); + return matchEverything; } DataDescriptorMatcher DataSpecUtils::dataDescriptorMatcherFrom(ConcreteDataTypeMatcher const& dataType) @@ -526,10 +526,10 @@ DataDescriptorMatcher DataSpecUtils::dataDescriptorMatcherFrom(ConcreteDataTypeM DescriptionValueMatcher{dataType.description.as()}, std::make_unique(DataDescriptorMatcher::Op::Just, StartTimeValueMatcher{ContextRef{0}})); - return std::move(DataDescriptorMatcher( + return DataDescriptorMatcher( DataDescriptorMatcher::Op::And, OriginValueMatcher{dataType.origin.as()}, - std::move(timeDescriptionMatcher))); + std::move(timeDescriptionMatcher)); } DataDescriptorMatcher DataSpecUtils::dataDescriptorMatcherFrom(header::DataOrigin const& origin) @@ -547,7 +547,7 @@ DataDescriptorMatcher DataSpecUtils::dataDescriptorMatcherFrom(header::DataOrigi SubSpecificationTypeValueMatcher{ContextRef{2}}, std::make_unique(DataDescriptorMatcher::Op::Just, StartTimeValueMatcher{ContextRef{0}})))}; - return std::move(matchOnlyOrigin); + return matchOnlyOrigin; } DataDescriptorMatcher DataSpecUtils::dataDescriptorMatcherFrom(header::DataDescription const& description) @@ -565,7 +565,7 @@ DataDescriptorMatcher DataSpecUtils::dataDescriptorMatcherFrom(header::DataDescr SubSpecificationTypeValueMatcher{ContextRef{2}}, std::make_unique(DataDescriptorMatcher::Op::Just, StartTimeValueMatcher{ContextRef{0}})))}; - return std::move(matchOnlyOrigin); + return matchOnlyOrigin; } std::optional DataSpecUtils::optionalConcreteDataMatcherFrom(data_matcher::DataDescriptorMatcher const& matcher) diff --git a/Framework/Core/src/DeviceSpecHelpers.cxx b/Framework/Core/src/DeviceSpecHelpers.cxx index c54dd6cec7a99..6c2063cb59311 100644 --- a/Framework/Core/src/DeviceSpecHelpers.cxx +++ b/Framework/Core/src/DeviceSpecHelpers.cxx @@ -48,6 +48,8 @@ #include #include +#include + #pragma GCC diagnostic push #pragma GCC diagnostic ignored "-Wpedantic" @@ -1072,6 +1074,8 @@ void DeviceSpecHelpers::dataProcessorSpecs2DeviceSpecs(const WorkflowSpec& workf WorkflowHelpers::constructGraph(workflow, logicalEdges, outputs, availableForwardsInfo); + WorkflowHelpers::validateEdges(workflow, logicalEdges, outputs); + // We need to instanciate one device per (me, timeIndex) in the // DeviceConnectionEdge. For each device we need one new binding // server per (me, other) -> port Moreover for each (me, other, diff --git a/Framework/Core/src/TopologyPolicy.cxx b/Framework/Core/src/TopologyPolicy.cxx index e0f7a208ed4e5..fb96eff5af2fc 100644 --- a/Framework/Core/src/TopologyPolicy.cxx +++ b/Framework/Core/src/TopologyPolicy.cxx @@ -82,33 +82,41 @@ bool expendableDataDeps(DataProcessorSpec const& a, DataProcessorSpec const& b) // If we are here we do not have any data dependency, // however we strill consider a dependent on b if // a has the "expendable" label and b does not. - bool isBExpendable = false; - bool isAExpendable = false; - for (auto const& label : b.labels) { + auto checkExpendable = [](DataProcessorLabel const& label) { if (label.value == "expendable") { - isBExpendable = true; - break; + return true; } - } - for (auto const& label : a.labels) { - if (label.value == "expendable") { - isAExpendable = true; - break; + return false; + }; + // A task marked as expendable or resilient can be put after an expendable task + auto checkResilient = [](DataProcessorLabel const& label) { + if (label.value == "resilient") { + return true; } - } - // If none is expendable. We simply return false. + return false; + }; + bool isBExpendable = std::find_if(b.labels.begin(), b.labels.end(), checkExpendable) != b.labels.end(); + bool isAExpendable = std::find_if(a.labels.begin(), a.labels.end(), checkExpendable) != a.labels.end(); + bool bResilient = std::find_if(b.labels.begin(), b.labels.end(), checkResilient) != b.labels.end(); + + // If none is expendable. We simply return false and sort as usual. if (!isAExpendable && !isBExpendable) { - LOGP(debug, "Neither {} nor {} are expendable. No dependency.", a.name, b.name); + LOGP(debug, "Neither {} nor {} are expendable. No dependency beyond data deps.", a.name, b.name); return false; } - // If both are expendable. We return false. + // If both are expendable. We return false and sort as usual. if (isAExpendable && isBExpendable) { LOGP(debug, "Both {} and {} are expendable. No dependency.", a.name, b.name); return false; } + // If b is expendable but b is resilient, we can keep the same order. + if (isAExpendable && bResilient) { + LOGP(debug, "{} is expendable but b is resilient, no need to add an unneeded dependency", a.name, a.name, b.name); + return false; + } // If a is expendable we consider it as if there was a dependency from a to b, - // but we still need to check if there is not one already from b to a. + // however we still need to check if there is not one already from b to a. if (isAExpendable) { LOGP(debug, "{} is expendable. Checking if there is a dependency from {} to {}.", a.name, b.name, a.name); return !dataDeps(b, a); diff --git a/Framework/Core/src/WorkflowHelpers.cxx b/Framework/Core/src/WorkflowHelpers.cxx index cd9c3913de228..2484f021ff0e2 100644 --- a/Framework/Core/src/WorkflowHelpers.cxx +++ b/Framework/Core/src/WorkflowHelpers.cxx @@ -11,7 +11,6 @@ #include "WorkflowHelpers.h" #include "Framework/AlgorithmSpec.h" #include "Framework/AODReaderHelpers.h" -#include "Framework/ChannelMatching.h" #include "Framework/ConfigParamsHelper.h" #include "Framework/CommonDataProcessors.h" #include "Framework/ConfigContext.h" @@ -21,13 +20,10 @@ #include "Framework/ControlService.h" #include "Framework/RawDeviceService.h" #include "Framework/StringHelpers.h" -#include "Framework/CommonMessageBackends.h" #include "Framework/ChannelSpecHelpers.h" -#include "Framework/ExternalFairMQDeviceProxy.h" #include "Framework/Plugins.h" #include "Framework/DataTakingContext.h" #include "Framework/DefaultsHelpers.h" -#include "ArrowSupport.h" #include "Headers/DataHeader.h" #include @@ -44,6 +40,10 @@ namespace o2::framework { +static constexpr std::array AODOrigins{header::DataOrigin{"AOD"}, header::DataOrigin{"AOD1"}, header::DataOrigin{"AOD2"}}; +static constexpr std::array extendedAODOrigins{header::DataOrigin{"AOD"}, header::DataOrigin{"AOD1"}, header::DataOrigin{"AOD2"}, header::DataOrigin{"DYN"}, header::DataOrigin{"AMD"}}; +static constexpr std::array writableAODOrigins{header::DataOrigin{"AOD"}, header::DataOrigin{"AOD1"}, header::DataOrigin{"AOD2"}, header::DataOrigin{"DYN"}}; + std::ostream& operator<<(std::ostream& out, TopoIndexInfo const& info) { out << "(" << info.index << ", " << info.layer << ")"; @@ -61,25 +61,25 @@ std::vector using EdgeIndex = int; // Create the index which will be returned. std::vector index(nodeCount); - for (auto wi = 0; wi < nodeCount; ++wi) { + for (auto wi = 0; static_cast(wi) < nodeCount; ++wi) { index[wi] = {wi, 0}; } std::vector remainingEdgesIndex(edgesCount); - for (EdgeIndex ei = 0; ei < edgesCount; ++ei) { + for (EdgeIndex ei = 0; static_cast(ei) < edgesCount; ++ei) { remainingEdgesIndex[ei] = ei; } // Create a vector where at each position we have true // if the vector has dependencies, false otherwise std::vector nodeDeps(nodeCount, false); - for (EdgeIndex ei = 0; ei < edgesCount; ++ei) { + for (EdgeIndex ei = 0; static_cast(ei) < edgesCount; ++ei) { nodeDeps[*(edgeOut + ei * stride)] = true; } // We start with all those which do not have any dependencies // They are layer 0. std::list L; - for (auto ii = 0; ii < index.size(); ++ii) { + for (auto ii = 0; static_cast(ii) < index.size(); ++ii) { if (nodeDeps[ii] == false) { L.push_back({ii, 0}); } @@ -201,7 +201,7 @@ void WorkflowHelpers::addMissingOutputsToBuilder(std::vector const& r if (j == publisher.inputs.end()) { publisher.inputs.push_back(spec); } - if (DataSpecUtils::partialMatch(spec, header::DataOrigin{"AOD"})) { + if (DataSpecUtils::partialMatch(spec, AODOrigins)) { DataSpecUtils::updateInputList(requestedAODs, std::move(spec)); } else if (DataSpecUtils::partialMatch(spec, header::DataOrigin{"DYN"})) { DataSpecUtils::updateInputList(requestedDYNs, std::move(spec)); @@ -424,7 +424,7 @@ void WorkflowHelpers::injectServiceDevices(WorkflowSpec& workflow, ConfigContext case Lifetime::Optional: break; } - if (DataSpecUtils::partialMatch(input, header::DataOrigin{"AOD"})) { + if (DataSpecUtils::partialMatch(input, AODOrigins)) { DataSpecUtils::updateInputList(requestedAODs, InputSpec{input}); } if (DataSpecUtils::partialMatch(input, header::DataOrigin{"DYN"})) { @@ -438,7 +438,7 @@ void WorkflowHelpers::injectServiceDevices(WorkflowSpec& workflow, ConfigContext std::stable_sort(timer.outputs.begin(), timer.outputs.end(), [](OutputSpec const& a, OutputSpec const& b) { return *DataSpecUtils::getOptionalSubSpec(a) < *DataSpecUtils::getOptionalSubSpec(b); }); for (auto& output : processor.outputs) { - if (DataSpecUtils::partialMatch(output, header::DataOrigin{"AOD"})) { + if (DataSpecUtils::partialMatch(output, AODOrigins)) { providedAODs.emplace_back(output); } else if (DataSpecUtils::partialMatch(output, header::DataOrigin{"DYN"})) { providedDYNs.emplace_back(output); @@ -670,13 +670,8 @@ void WorkflowHelpers::injectServiceDevices(WorkflowSpec& workflow, ConfigContext // ATTENTION: if there are dangling outputs the getGlobalAODSink // has to be created in any case! std::vector outputsInputsAOD; - auto isAOD = [](InputSpec const& spec) { - return (DataSpecUtils::partialMatch(spec, header::DataOrigin("AOD")) || - DataSpecUtils::partialMatch(spec, header::DataOrigin("DYN")) || - DataSpecUtils::partialMatch(spec, header::DataOrigin("AMD"))); - }; for (auto ii = 0u; ii < outputsInputs.size(); ii++) { - if (isAOD(outputsInputs[ii])) { + if (DataSpecUtils::partialMatch(outputsInputs[ii], extendedAODOrigins)) { auto ds = dod->getDataOutputDescriptors(outputsInputs[ii]); if (ds.size() > 0 || isDangling[ii]) { outputsInputsAOD.emplace_back(outputsInputs[ii]); @@ -714,7 +709,7 @@ void WorkflowHelpers::injectServiceDevices(WorkflowSpec& workflow, ConfigContext continue; } // AODs are skipped in any case. - if (isAOD(outputsInputs[ii])) { + if (DataSpecUtils::partialMatch(outputsInputs[ii], extendedAODOrigins)) { continue; } redirectedOutputsInputs.emplace_back(outputsInputs[ii]); @@ -1118,26 +1113,21 @@ std::shared_ptr WorkflowHelpers::getDataOutputDirector(Confi } } // parse the keepString - auto isAOD = [](InputSpec const& spec) { return DataSpecUtils::partialMatch(spec, header::DataOrigin("AOD")); }; if (options.isSet("aod-writer-keep")) { auto keepString = options.get("aod-writer-keep"); if (!keepString.empty()) { - dod->reset(); std::string d("dangling"); if (d.find(keepString) == 0) { - // use the dangling outputs std::vector danglingOutputs; for (auto ii = 0u; ii < OutputsInputs.size(); ii++) { - if (isAOD(OutputsInputs[ii]) && isDangling[ii]) { + if (DataSpecUtils::partialMatch(OutputsInputs[ii], writableAODOrigins) && isDangling[ii]) { danglingOutputs.emplace_back(OutputsInputs[ii]); } } dod->readSpecs(danglingOutputs); - } else { - // use the keep string dod->readString(keepString); } @@ -1232,5 +1222,59 @@ std::vector WorkflowHelpers::computeDanglingOutputs(WorkflowSpec cons return results; } -#pragma diagnostic pop +bool validateLifetime(std::ostream& errors, DataProcessorSpec const& producer, OutputSpec const& output, DataProcessorSpec const& consumer, InputSpec const& input) +{ + if (input.lifetime == Lifetime::Timeframe && output.lifetime == Lifetime::Sporadic) { + errors << fmt::format("Input {} of {} has lifetime Timeframe, but output {} of {} has lifetime Sporadic\n", + DataSpecUtils::describe(input).c_str(), consumer.name, + DataSpecUtils::describe(output).c_str(), producer.name); + return false; + } + return true; +} + +bool validateExpendable(std::ostream& errors, DataProcessorSpec const& producer, OutputSpec const& output, DataProcessorSpec const& consumer, InputSpec const& input) +{ + auto isExpendable = [](DataProcessorLabel const& label) { + return label.value == "expendable"; + }; + auto isResilient = [](DataProcessorLabel const& label) { + return label.value == "expendable" || label.value == "resilient"; + }; + bool producerExpendable = std::find_if(producer.labels.begin(), producer.labels.end(), isExpendable) != producer.labels.end(); + bool consumerCritical = std::find_if(consumer.labels.begin(), consumer.labels.end(), isResilient) == consumer.labels.end(); + if (producerExpendable && consumerCritical) { + errors << fmt::format("Critical consumer {} depends on expendable producer {}\n", + consumer.name, + producer.name); + return false; + } + return true; +} + +using Validator = std::function; +void WorkflowHelpers::validateEdges(WorkflowSpec const& workflow, + std::vector const& edges, + std::vector const& outputs) +{ + std::vector defaultValidators = {validateExpendable, validateLifetime}; + std::stringstream errors; + // Iterate over all the edges. + // Get the input lifetime and the output lifetime. + // Output lifetime must be Timeframe if the input lifetime is Timeframe. + bool hasErrors = false; + for (auto& edge : edges) { + DataProcessorSpec const& producer = workflow[edge.producer]; + DataProcessorSpec const& consumer = workflow[edge.consumer]; + OutputSpec const& output = outputs[edge.outputGlobalIndex]; + InputSpec const& input = consumer.inputs[edge.consumerInputIndex]; + for (auto& validator : defaultValidators) { + hasErrors |= !validator(errors, producer, output, consumer, input); + } + } + if (hasErrors) { + throw std::runtime_error(errors.str()); + } +} + } // namespace o2::framework diff --git a/Framework/Core/src/WorkflowHelpers.h b/Framework/Core/src/WorkflowHelpers.h index a167e99d4768c..64bd593b68697 100644 --- a/Framework/Core/src/WorkflowHelpers.h +++ b/Framework/Core/src/WorkflowHelpers.h @@ -13,10 +13,8 @@ #include "Framework/InputSpec.h" #include "Framework/OutputSpec.h" -#include "Framework/ForwardRoute.h" #include "Framework/WorkflowSpec.h" #include "Framework/DataOutputDirector.h" -#include "Framework/DataProcessorInfo.h" #include #include @@ -227,6 +225,14 @@ struct WorkflowHelpers { /// returns only dangling outputs static std::vector computeDanglingOutputs(WorkflowSpec const& workflow); + + /// Validate that the nodes at the ends of the edges of the graph + /// are actually compatible with each other. + /// For example we should make sure that Lifetime::Timeframe inputs of + /// one node is not connected to an Output of Lifetime::Sporadic of another node. + static void validateEdges(WorkflowSpec const& workflow, + std::vector const& edges, + std::vector const& outputs); }; } // namespace o2::framework diff --git a/Framework/Core/src/runDataProcessing.cxx b/Framework/Core/src/runDataProcessing.cxx index 86cc15299e8f3..724383ce501e6 100644 --- a/Framework/Core/src/runDataProcessing.cxx +++ b/Framework/Core/src/runDataProcessing.cxx @@ -2565,6 +2565,30 @@ void apply_permutation( } } +// Check if the workflow is resiliant to failures +void checkNonResiliency(std::vector const& specs, + std::vector> const& edges) +{ + auto checkExpendable = [](DataProcessorLabel const& label) { + return label.value == "expendable"; + }; + auto checkResilient = [](DataProcessorLabel const& label) { + return label.value == "resilient" || label.value == "expendable"; + }; + + for (auto& edge : edges) { + auto& src = specs[edge.first]; + auto& dst = specs[edge.second]; + if (std::none_of(src.labels.begin(), src.labels.end(), checkExpendable)) { + continue; + } + if (std::any_of(dst.labels.begin(), dst.labels.end(), checkResilient)) { + continue; + } + throw std::runtime_error("Workflow is not resiliant to failures. Processor " + dst.name + " gets inputs from expendable devices, but is not marked as expendable or resilient itself."); + } +} + std::string debugTopoInfo(std::vector const& specs, std::vector const& infos, std::vector> const& edges) @@ -2590,6 +2614,11 @@ std::string debugTopoInfo(std::vector const& specs, for (auto& d : specs) { out << "- " << d.name << std::endl; } + out << "digraph G {\n"; + for (auto& e : edges) { + out << fmt::format(" \"{}\" -> \"{}\"\n", specs[e.first].name, specs[e.second].name); + } + out << "}\n"; return out.str(); } @@ -2828,6 +2857,8 @@ int doMain(int argc, char** argv, o2::framework::WorkflowSpec const& workflow, auto topoInfos = WorkflowHelpers::topologicalSort(physicalWorkflow.size(), &edges[0].first, &edges[0].second, sizeof(std::pair), edges.size()); if (topoInfos.size() != physicalWorkflow.size()) { + // Check missing resilincy of one of the tasks + checkNonResiliency(physicalWorkflow, edges); throw std::runtime_error("Unable to do topological sort of the resulting workflow. Do you have loops?\n" + debugTopoInfo(physicalWorkflow, topoInfos, edges)); } // Sort by layer and then by name, to ensure stability. diff --git a/Framework/Core/test/test_FrameworkDataFlowToO2Control.cxx b/Framework/Core/test/test_FrameworkDataFlowToO2Control.cxx index 4ef883f3b32de..7b0ffa462a9db 100644 --- a/Framework/Core/test/test_FrameworkDataFlowToO2Control.cxx +++ b/Framework/Core/test/test_FrameworkDataFlowToO2Control.cxx @@ -31,32 +31,26 @@ namespace { WorkflowSpec defineDataProcessing() { - return {{"A", // - Inputs{}, // - Outputs{OutputSpec{"TST", "A1"}, OutputSpec{"TST", "A2"}}, // A1 will be consumed twice, A2 is dangling - AlgorithmSpec{}, // - {ConfigParamSpec{"channel-config", VariantType::String, // raw input channel - "name=into_dpl,type=pull,method=connect,address=ipc:///tmp/pipe-into-dpl,transport=shmem,rateLogging=10,rcvBufSize=789", - {"Out-of-band channel config"}}}}, - {"B", // producer, no inputs - Inputs{}, - Outputs{OutputSpec{"TST", "B1"}}, + return {{.name = "A", // + .outputs = Outputs{OutputSpec{"TST", "A1"}, OutputSpec{"TST", "A2"}}, // A1 will be consumed twice, A2 is dangling + .algorithm = AlgorithmSpec{}, // + .options = {ConfigParamSpec{"channel-config", VariantType::String, // raw input channel + "name=into_dpl,type=pull,method=connect,address=ipc:///tmp/pipe-into-dpl,transport=shmem,rateLogging=10,rcvBufSize=789", + {"Out-of-band channel config"}}}}, + {.name = "B", // producer, no inputs + .outputs = Outputs{OutputSpec{"TST", "B1"}}, .metadata = {{ecs::cpuKillThreshold, "3.0"}}}, - {"C", // first consumer of A1, consumer of B1 - {InputSpec{"y", "TST", "A1"}, InputSpec{"y", "TST", "B1"}}, - Outputs{}, + {.name = "C", // first consumer of A1, consumer of B1 + .inputs = {InputSpec{"y", "TST", "A1"}, InputSpec{"y", "TST", "B1"}}, .metadata = {{ecs::privateMemoryKillThresholdMB, "5000"}}}, - {"D", // second consumer of A1 - Inputs{ - InputSpec{"x", "TST", "A1"}}, - Outputs{}, - AlgorithmSpec{}, - {ConfigParamSpec{"a-param", VariantType::Int, 1, {"A parameter which should not be escaped"}}, - ConfigParamSpec{"b-param", VariantType::String, "", {"a parameter which will be escaped"}}, - ConfigParamSpec{"c-param", VariantType::String, "foo;bar", {"another parameter which will be escaped"}}, - ConfigParamSpec{"channel-config", VariantType::String, // raw output channel - "name=outta_dpl,type=push,method=bind,address=ipc:///tmp/pipe-outta-dpl,transport=shmem,rateLogging=10", - {"Out-of-band channel config"}}}}}; + {.name = "D", // second consumer of A1 + .inputs = Inputs{InputSpec{"x", "TST", "A1"}}, + .options = {ConfigParamSpec{"a-param", VariantType::Int, 1, {"A parameter which should not be escaped"}}, + ConfigParamSpec{"b-param", VariantType::String, "", {"a parameter which will be escaped"}}, + ConfigParamSpec{"c-param", VariantType::String, "foo;bar", {"another parameter which will be escaped"}}, + ConfigParamSpec{"channel-config", VariantType::String, // raw output channel + "name=outta_dpl,type=push,method=bind,address=ipc:///tmp/pipe-outta-dpl,transport=shmem,rateLogging=10", + {"Out-of-band channel config"}}}}}; } char* strdiffchr(const char* s1, const char* s2) diff --git a/GPU/GPUTracking/DataTypes/GPUTRDInterfaceO2Track.h b/GPU/GPUTracking/DataTypes/GPUTRDInterfaceO2Track.h index 19cad3649fa42..fb8f6d1b47d52 100644 --- a/GPU/GPUTracking/DataTypes/GPUTRDInterfaceO2Track.h +++ b/GPU/GPUTracking/DataTypes/GPUTRDInterfaceO2Track.h @@ -47,10 +47,10 @@ template <> class trackInterface : public o2::track::TrackParCov { public: - GPUdDefault() trackInterface() = default; - trackInterface(const o2::track::TrackParCov& param) = delete; - GPUd() trackInterface(const o2::dataformats::TrackTPCITS& trkItsTpc) : o2::track::TrackParCov(trkItsTpc.getParamOut()) {} - GPUd() trackInterface(const o2::tpc::TrackTPC& trkTpc) : o2::track::TrackParCov(trkTpc.getParamOut()) {} + GPUdDefault() trackInterface() = default; + trackInterface(const o2::track::TrackParCov& param) = delete; + GPUd() trackInterface(const o2::dataformats::TrackTPCITS& trkItsTpc) : o2::track::TrackParCov(trkItsTpc.getParamOut()) {} + GPUd() trackInterface(const o2::tpc::TrackTPC& trkTpc) : o2::track::TrackParCov(trkTpc.getParamOut()) {} GPUd() void set(float x, float alpha, const float* param, const float* cov) { @@ -63,8 +63,8 @@ class trackInterface : public o2::track::TrackParCov setCov(cov[i], i); } } - GPUd() trackInterface(const GPUTPCGMMergedTrack& trk); - GPUd() trackInterface(const gputpcgmmergertypes::GPUTPCOuterParam& param); + GPUd() trackInterface(const GPUTPCGMMergedTrack& trk); + GPUd() trackInterface(const gputpcgmmergertypes::GPUTPCOuterParam& param); GPUd() void updateCovZ2(float addZerror) { updateCov(addZerror, o2::track::CovLabels::kSigZ2); } GPUd() o2::track::TrackLTIntegral& getLTIntegralOut() { return mLTOut; } GPUd() const o2::track::TrackLTIntegral& getLTIntegralOut() const { return mLTOut; } diff --git a/GPU/GPUTracking/TRDTracking/GPUTRDInterfaces.h b/GPU/GPUTracking/TRDTracking/GPUTRDInterfaces.h index ab8db40e90e2e..6359f880e2860 100644 --- a/GPU/GPUTracking/TRDTracking/GPUTRDInterfaces.h +++ b/GPU/GPUTracking/TRDTracking/GPUTRDInterfaces.h @@ -101,9 +101,9 @@ class propagatorInterface : public AliTrackerBase public: typedef void propagatorParam; - propagatorInterface(const propagatorParam* = nullptr) : AliTrackerBase(), mParam(nullptr){}; - propagatorInterface(const propagatorInterface&) CON_DELETE; - propagatorInterface& operator=(const propagatorInterface&) CON_DELETE; + propagatorInterface(const propagatorParam* = nullptr) : AliTrackerBase(), mParam(nullptr){}; + propagatorInterface(const propagatorInterface&) CON_DELETE; + propagatorInterface& operator=(const propagatorInterface&) CON_DELETE; bool propagateToX(float x, float maxSnp, float maxStep) { return PropagateTrackToBxByBz(mParam, x, 0.13957, maxStep, false, maxSnp); } int getPropagatedYZ(float x, float& projY, float& projZ) @@ -147,9 +147,9 @@ class propagatorInterface { public: typedef o2::base::Propagator propagatorParam; - GPUd() propagatorInterface(const propagatorParam* prop) : mProp(prop){}; - GPUd() propagatorInterface(const propagatorInterface&) = delete; - GPUd() propagatorInterface& operator=(const propagatorInterface&) = delete; + GPUd() propagatorInterface(const propagatorParam* prop) : mProp(prop){}; + GPUd() propagatorInterface(const propagatorInterface&) = delete; + GPUd() propagatorInterface& operator=(const propagatorInterface&) = delete; GPUdi() bool propagateToX(float x, float maxSnp, float maxStep) { return mProp->PropagateToXBxByBz(*mParam, x, maxSnp, maxStep); } GPUdi() int getPropagatedYZ(float x, float& projY, float& projZ) { return static_cast(mParam->getYZAt(x, mProp->getNominalBz(), projY, projZ)); } @@ -206,10 +206,10 @@ template <> class trackInterface : public GPUTPCGMTrackParam { public: - GPUdDefault() trackInterface() CON_DEFAULT; - GPUd() trackInterface(const GPUTPCGMTrackParam& param) CON_DELETE; - GPUd() trackInterface(const GPUTPCGMMergedTrack& trk) : GPUTPCGMTrackParam(trk.GetParam()), mAlpha(trk.GetAlpha()) {} - GPUd() trackInterface(const gputpcgmmergertypes::GPUTPCOuterParam& param) : GPUTPCGMTrackParam(), mAlpha(param.alpha) + GPUdDefault() trackInterface() CON_DEFAULT; + GPUd() trackInterface(const GPUTPCGMTrackParam& param) CON_DELETE; + GPUd() trackInterface(const GPUTPCGMMergedTrack& trk) : GPUTPCGMTrackParam(trk.GetParam()), mAlpha(trk.GetAlpha()) {} + GPUd() trackInterface(const gputpcgmmergertypes::GPUTPCOuterParam& param) : GPUTPCGMTrackParam(), mAlpha(param.alpha) { SetX(param.X); for (int i = 0; i < 5; i++) { @@ -220,11 +220,11 @@ class trackInterface : public GPUTPCGMTrackParam } }; #ifdef GPUCA_NOCOMPAT - GPUdDefault() trackInterface(const trackInterface& param) = default; - GPUdDefault() trackInterface& operator=(const trackInterface& param) = default; + GPUdDefault() trackInterface(const trackInterface& param) = default; + GPUdDefault() trackInterface& operator=(const trackInterface& param) = default; #endif #ifdef GPUCA_ALIROOT_LIB - trackInterface(const AliHLTExternalTrackParam& param) : GPUTPCGMTrackParam(), mAlpha(param.fAlpha) + trackInterface(const AliHLTExternalTrackParam& param) : GPUTPCGMTrackParam(), mAlpha(param.fAlpha) { SetX(param.fX); SetPar(0, param.fY); @@ -238,7 +238,7 @@ class trackInterface : public GPUTPCGMTrackParam }; #endif #if defined(GPUCA_HAVE_O2HEADERS) - GPUd() trackInterface(const o2::dataformats::TrackTPCITS& param) : GPUTPCGMTrackParam(), mAlpha(param.getParamOut().getAlpha()) + GPUd() trackInterface(const o2::dataformats::TrackTPCITS& param) : GPUTPCGMTrackParam(), mAlpha(param.getParamOut().getAlpha()) { SetX(param.getParamOut().getX()); SetPar(0, param.getParamOut().getY()); @@ -250,7 +250,7 @@ class trackInterface : public GPUTPCGMTrackParam SetCov(i, param.getParamOut().getCov()[i]); } } - GPUd() trackInterface(const o2::tpc::TrackTPC& param) : GPUTPCGMTrackParam(), mAlpha(param.getParamOut().getAlpha()) + GPUd() trackInterface(const o2::tpc::TrackTPC& param) : GPUTPCGMTrackParam(), mAlpha(param.getParamOut().getAlpha()) { SetX(param.getParamOut().getX()); SetPar(0, param.getParamOut().getY()); @@ -307,7 +307,7 @@ class propagatorInterface : public GPUTPCGMPropagator { public: typedef GPUTPCGMPolynomialField propagatorParam; - GPUd() propagatorInterface(const propagatorParam* pField) : GPUTPCGMPropagator(), mTrack(nullptr) + GPUd() propagatorInterface(const propagatorParam* pField) : GPUTPCGMPropagator(), mTrack(nullptr) { this->SetMaterialTPC(); this->SetPolynomialField(pField); @@ -316,8 +316,8 @@ class propagatorInterface : public GPUTPCGMPropagator this->SetFitInProjections(0); this->SelectFieldRegion(GPUTPCGMPropagator::TRD); }; - propagatorInterface(const propagatorInterface&) CON_DELETE; - propagatorInterface& operator=(const propagatorInterface&) CON_DELETE; + propagatorInterface(const propagatorInterface&) CON_DELETE; + propagatorInterface& operator=(const propagatorInterface&) CON_DELETE; GPUd() void setTrack(trackInterface* trk) { SetTrack(trk, trk->getAlpha()); diff --git a/GPU/TPCFastTransformation/CorrectionMapsHelper.h b/GPU/TPCFastTransformation/CorrectionMapsHelper.h index 75d194bef63cb..970f9b4c88a26 100644 --- a/GPU/TPCFastTransformation/CorrectionMapsHelper.h +++ b/GPU/TPCFastTransformation/CorrectionMapsHelper.h @@ -99,7 +99,7 @@ class CorrectionMapsHelper { if (mMeanLumi < 0.f || mInstLumi < 0.f) { mLumiScale = -1.f; - } else if (mLumiScaleMode == 1) { + } else if ((mLumiScaleMode == 1) || (mLumiScaleMode == 2)) { mLumiScale = mMeanLumiRef ? (mInstLumi - mMeanLumi) / mMeanLumiRef : 0.f; LOGP(debug, "mInstLumi: {} mMeanLumi: {} mMeanLumiRef: {}", mInstLumi, mMeanLumi, mMeanLumiRef); } else { diff --git a/GPU/TPCFastTransformation/TPCFastTransform.h b/GPU/TPCFastTransformation/TPCFastTransform.h index 956ae13646ff6..1936ae4f19a99 100644 --- a/GPU/TPCFastTransformation/TPCFastTransform.h +++ b/GPU/TPCFastTransformation/TPCFastTransform.h @@ -442,7 +442,7 @@ GPUdi() void TPCFastTransform::TransformInternal(int slice, int row, float& u, f { if (mApplyCorrection) { float dx = 0.f, du = 0.f, dv = 0.f; - if ((scale >= 0.f) || (scaleMode == 1)) { + if ((scale >= 0.f) || (scaleMode == 1) || (scaleMode == 2)) { #ifndef GPUCA_GPUCODE if (mCorrectionSlow) { float ly, lz; @@ -470,7 +470,7 @@ GPUdi() void TPCFastTransform::TransformInternal(int slice, int row, float& u, f dx = (dx - dxRef) * scale + dxRef; du = (du - duRef) * scale + duRef; dv = (dv - dvRef) * scale + dvRef; - } else if ((scale != 0.f) && (scaleMode == 1)) { + } else if ((scale != 0.f) && ((scaleMode == 1) || (scaleMode == 2))) { float dxRef, duRef, dvRef; ref->mCorrection.getCorrection(slice, row, u, v, dxRef, duRef, dvRef); dx = dxRef * scale + dx; @@ -494,11 +494,11 @@ GPUdi() void TPCFastTransform::TransformInternal(int slice, int row, float& u, f getGeometry().convUVtoLocal(slice, uCorr, vCorr, lyT, lzT); float invYZtoX; - InverseTransformYZtoX(slice, row, ly, lz, invYZtoX); + InverseTransformYZtoX(slice, row, lyT, lzT, invYZtoX, ref, scale, scaleMode); float YZtoNominalY; float YZtoNominalZ; - InverseTransformYZtoNominalYZ(slice, row, ly, lz, YZtoNominalY, YZtoNominalZ); + InverseTransformYZtoNominalYZ(slice, row, lyT, lzT, YZtoNominalY, YZtoNominalZ, ref, scale, scaleMode); float dxRef, duRef, dvRef; ref->mCorrection.getCorrection(slice, row, u, v, dxRef, duRef, dvRef); @@ -745,17 +745,17 @@ GPUdi() void TPCFastTransform::InverseTransformYZtoX(int slice, int row, float y /// Transformation y,z -> x float u = 0, v = 0; getGeometry().convLocalToUV(slice, y, z, u, v); - if (scale >= 0.f) { + if ((scale >= 0.f) || (scaleMode == 1) || (scaleMode == 2)) { mCorrection.getCorrectionInvCorrectedX(slice, row, u, v, x); - if (ref && scale > 0.f) { // scaling was requested - if (scaleMode == 0) { + if (ref) { // scaling was requested + if (scaleMode == 0 && scale > 0.f) { float xr; ref->mCorrection.getCorrectionInvCorrectedX(slice, row, u, v, xr); x = (x - xr) * scale + xr; - } else if (scaleMode == 1) { + } else if ((scale != 0) && ((scaleMode == 1) || (scaleMode == 2))) { float xr; ref->mCorrection.getCorrectionInvCorrectedX(slice, row, u, v, xr); - x = xr * scale + x; + x = (xr - getGeometry().getRowInfo(row).x) * scale + x; // xr=mGeo.getRowInfo(row).x + dx; } } } else { @@ -780,19 +780,19 @@ GPUdi() void TPCFastTransform::InverseTransformYZtoNominalYZ(int slice, int row, /// Transformation y,z -> x float u = 0, v = 0, un = 0, vn = 0; getGeometry().convLocalToUV(slice, y, z, u, v); - if (scale >= 0.f) { + if ((scale >= 0.f) || (scaleMode == 1) || (scaleMode == 2)) { mCorrection.getCorrectionInvUV(slice, row, u, v, un, vn); - if (ref && scale > 0.f) { // scaling was requested - if (scaleMode == 0) { + if (ref) { // scaling was requested + if (scaleMode == 0 && scale > 0.f) { float unr = 0, vnr = 0; ref->mCorrection.getCorrectionInvUV(slice, row, u, v, unr, vnr); un = (un - unr) * scale + unr; vn = (vn - vnr) * scale + vnr; - } else if (scaleMode == 1) { + } else if ((scale != 0) && ((scaleMode == 1) || (scaleMode == 2))) { float unr = 0, vnr = 0; ref->mCorrection.getCorrectionInvUV(slice, row, u, v, unr, vnr); - un = unr * scale + un; - vn = vnr * scale + vn; + un = (unr - u) * scale + un; // unr = u - duv[0]; + vn = (vnr - v) * scale + vn; } } } else { diff --git a/Steer/DigitizerWorkflow/src/SimpleDigitizerWorkflow.cxx b/Steer/DigitizerWorkflow/src/SimpleDigitizerWorkflow.cxx index 79555249dd64c..f1e3bdff9a0c6 100644 --- a/Steer/DigitizerWorkflow/src/SimpleDigitizerWorkflow.cxx +++ b/Steer/DigitizerWorkflow/src/SimpleDigitizerWorkflow.cxx @@ -169,6 +169,7 @@ void customize(std::vector& workflowOptions) // Option to write TPC digits internaly, without forwarding to a special writer instance. // This is useful in GRID productions with small available memory. workflowOptions.push_back(ConfigParamSpec{"tpc-chunked-writer", o2::framework::VariantType::Bool, false, {"Write independent TPC digit chunks as soon as they can be flushed."}}); + workflowOptions.push_back(ConfigParamSpec{"tpc-distortion-type", o2::framework::VariantType::Int, 0, {"Simulate distortions in the TPC (0=no distortions, 1=distortions without scaling, 2=distortions with CTP scaling)"}}); std::string simhelp("Comma separated list of simulation prefixes (for background, signal productions)"); workflowOptions.push_back( @@ -568,7 +569,8 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) detList.emplace_back(o2::detectors::DetID::TPC); auto internalwrite = configcontext.options().get("tpc-chunked-writer"); - WorkflowSpec tpcPipelines = o2::tpc::getTPCDigitizerSpec(lanes, tpcsectors, mctruth, internalwrite); + auto distortionType = configcontext.options().get("tpc-distortion-type"); + WorkflowSpec tpcPipelines = o2::tpc::getTPCDigitizerSpec(lanes, tpcsectors, mctruth, internalwrite, distortionType); specs.insert(specs.end(), tpcPipelines.begin(), tpcPipelines.end()); if (configcontext.options().get("tpc-reco-type").empty() == false) { diff --git a/Steer/DigitizerWorkflow/src/TPCDigitizerSpec.cxx b/Steer/DigitizerWorkflow/src/TPCDigitizerSpec.cxx index 1b576008fdcd6..59f8589a56449 100644 --- a/Steer/DigitizerWorkflow/src/TPCDigitizerSpec.cxx +++ b/Steer/DigitizerWorkflow/src/TPCDigitizerSpec.cxx @@ -44,12 +44,12 @@ #include "CommonDataFormat/RangeReference.h" #include "SimConfig/DigiParams.h" #include -#include "TH3.h" +#include "Framework/CCDBParamSpec.h" using namespace o2::framework; using SubSpecificationType = o2::framework::DataAllocator::SubSpecificationType; using DigiGroupRef = o2::dataformats::RangeReference; -using SC = o2::tpc::SpaceCharge; +using SC = o2::tpc::SpaceCharge; namespace o2 { @@ -108,7 +108,7 @@ using namespace o2::base; class TPCDPLDigitizerTask : public BaseDPLDigitizer { public: - TPCDPLDigitizerTask(bool internalwriter) : mInternalWriter(internalwriter), BaseDPLDigitizer(InitServices::FIELD | InitServices::GEOM) + TPCDPLDigitizerTask(bool internalwriter, int distortionType) : mInternalWriter(internalwriter), BaseDPLDigitizer(InitServices::FIELD | InitServices::GEOM), mDistortionType(distortionType) { } @@ -119,61 +119,15 @@ class TPCDPLDigitizerTask : public BaseDPLDigitizer mLaneId = ic.services().get().rank; mWithMCTruth = o2::conf::DigiParams::Instance().mctruth; - auto useDistortions = ic.options().get("distortionType"); auto triggeredMode = ic.options().get("TPCtriggered"); mUseCalibrationsFromCCDB = ic.options().get("TPCuseCCDB"); + mMeanLumiDistortions = ic.options().get("meanLumiDistortions"); + mMeanLumiDistortionsDerivative = ic.options().get("meanLumiDistortionsDerivative"); + LOG(info) << "TPC calibrations from CCDB: " << mUseCalibrationsFromCCDB; - if (useDistortions > 0) { - if (useDistortions == 1) { - LOG(info) << "Using realistic space-charge distortions."; - } else { - LOG(info) << "Using constant space-charge distortions."; - } - auto readSpaceChargeString = ic.options().get("readSpaceCharge"); - std::vector readSpaceCharge; - std::stringstream ssSpaceCharge(readSpaceChargeString); - while (ssSpaceCharge.good()) { - std::string substr; - getline(ssSpaceCharge, substr, ','); - readSpaceCharge.push_back(substr); - } - if (readSpaceCharge[0].size() != 0) { // use pre-calculated space-charge object - if (std::filesystem::exists(readSpaceCharge[0])) { - LOGP(info, "Reading space-charge object from file {}", readSpaceCharge[0].data()); - mDigitizer.setUseSCDistortions(readSpaceCharge[0]); - } else { - LOG(error) << "Space-charge object or file not found!"; - } - } else { // create new space-charge object either with empty TPC or an initial space-charge density provided by histogram - SCDistortionType distortionType = useDistortions == 2 ? SCDistortionType::SCDistortionsConstant : SCDistortionType::SCDistortionsRealistic; - auto inputHistoString = ic.options().get("initialSpaceChargeDensity"); - std::vector inputHisto; - std::stringstream ssHisto(inputHistoString); - while (ssHisto.good()) { - std::string substr; - getline(ssHisto, substr, ','); - inputHisto.push_back(substr); - } - std::unique_ptr hisSCDensity; - if (std::filesystem::exists(inputHisto[0])) { - auto fileSCInput = std::unique_ptr(TFile::Open(inputHisto[0].data())); - if (fileSCInput->FindKey(inputHisto[1].data())) { - hisSCDensity.reset((TH3*)fileSCInput->Get(inputHisto[1].data())); - hisSCDensity->SetDirectory(nullptr); - } - } - if (hisSCDensity.get() != nullptr) { - LOG(info) << "TPC: Providing initial space-charge density histogram: " << hisSCDensity->GetName(); - mDigitizer.setUseSCDistortions(distortionType, hisSCDensity.get()); - } else { - if (distortionType == SCDistortionType::SCDistortionsConstant) { - LOG(error) << "Input space-charge density histogram or file not found!"; - } - } - } - } mDigitizer.setContinuousReadout(!triggeredMode); + mDigitizer.setDistortionScaleType(mDistortionType); // we send the GRP data once if the corresponding output channel is available // and set the flag to false after @@ -236,6 +190,20 @@ class TPCDPLDigitizerTask : public BaseDPLDigitizer if (mTPCVDriftHelper.accountCCDBInputs(matcher, obj)) { return; } + if (matcher == ConcreteDataMatcher(o2::header::gDataOriginTPC, "TPCDIST", 0)) { + LOGP(info, "Updating distortion map"); + mDigitizer.setUseSCDistortions(static_cast(obj)); + if (mMeanLumiDistortions >= 0) { + mDigitizer.setMeanLumiDistortions(mMeanLumiDistortions); + } + } + if (matcher == ConcreteDataMatcher(o2::header::gDataOriginTPC, "TPCDISTDERIV", 0)) { + LOGP(info, "Updating reference distortion map"); + mDigitizer.setSCDistortionsDerivative(static_cast(obj)); + if (mMeanLumiDistortionsDerivative >= 0) { + mDigitizer.setMeanLumiDistortionsDerivative(mMeanLumiDistortionsDerivative); + } + } } void run(framework::ProcessingContext& pc) @@ -247,6 +215,14 @@ class TPCDPLDigitizerTask : public BaseDPLDigitizer cdb.setUseDefaults(!mUseCalibrationsFromCCDB); // whatever are global settings for CCDB usage, we have to extract the TPC vdrift from CCDB for anchored simulations mTPCVDriftHelper.extractCCDBInputs(pc); + if (mDistortionType) { + pc.inputs().get("tpcdistortions"); + if (mDistortionType == 2) { + pc.inputs().get("tpcdistortionsderiv"); + mDigitizer.setLumiScaleFactor(); + } + } + if (mTPCVDriftHelper.isUpdated()) { const auto& vd = mTPCVDriftHelper.getVDriftObject(); LOGP(info, "Updating TPC fast transform map with new VDrift factor of {} wrt reference {} and DriftTimeOffset correction {} wrt {} from source {}", @@ -498,9 +474,12 @@ class TPCDPLDigitizerTask : public BaseDPLDigitizer bool mWithMCTruth = true; bool mInternalWriter = false; bool mUseCalibrationsFromCCDB = false; + int mDistortionType = 0; + float mMeanLumiDistortions = -1; + float mMeanLumiDistortionsDerivative = -1; }; -o2::framework::DataProcessorSpec getTPCDigitizerSpec(int channel, bool writeGRP, bool mctruth, bool internalwriter) +o2::framework::DataProcessorSpec getTPCDigitizerSpec(int channel, bool writeGRP, bool mctruth, bool internalwriter, int distortionType) { // create the full data processor spec using // a name identifier @@ -530,23 +509,22 @@ o2::framework::DataProcessorSpec getTPCDigitizerSpec(int channel, bool writeGRP, id.str().c_str(), inputs, outputs, - AlgorithmSpec{adaptFromTask(internalwriter)}, + AlgorithmSpec{adaptFromTask(internalwriter, distortionType)}, Options{ - {"distortionType", VariantType::Int, 0, {"Distortion type to be used. 0 = no distortions (default), 1 = realistic distortions (not implemented yet), 2 = constant distortions"}}, - {"initialSpaceChargeDensity", VariantType::String, "", {"Path to root file containing TH3 with initial space-charge density and name of the TH3 (comma separated)"}}, - {"readSpaceCharge", VariantType::String, "", {"Path to root file containing pre-calculated space-charge object and name of the object (comma separated)"}}, {"TPCtriggered", VariantType::Bool, false, {"Impose triggered RO mode (default: continuous)"}}, {"TPCuseCCDB", VariantType::Bool, false, {"true: load calibrations from CCDB; false: use random calibratoins"}}, + {"meanLumiDistortions", VariantType::Float, -1.f, {"override lumi of distortion object if >=0"}}, + {"meanLumiDistortionsDerivative", VariantType::Float, -1.f, {"override lumi of derivative distortion object if >=0"}}, }}; } -o2::framework::WorkflowSpec getTPCDigitizerSpec(int nLanes, std::vector const& sectors, bool mctruth, bool internalwriter) +o2::framework::WorkflowSpec getTPCDigitizerSpec(int nLanes, std::vector const& sectors, bool mctruth, bool internalwriter, int distortionType) { // channel parameter is deprecated in the TPCDigitizer processor, all descendants // are initialized not to publish GRP mode, but the channel will be added to the first // processor after the pipelines have been created. The processor will decide upon // the index in the ParallelContext whether to publish - WorkflowSpec pipelineTemplate{getTPCDigitizerSpec(0, false, mctruth, internalwriter)}; + WorkflowSpec pipelineTemplate{getTPCDigitizerSpec(0, false, mctruth, internalwriter, distortionType)}; // override the predefined name, index will be added by parallelPipeline method pipelineTemplate[0].name = "TPCDigitizer"; WorkflowSpec pipelines = parallelPipeline( @@ -554,6 +532,13 @@ o2::framework::WorkflowSpec getTPCDigitizerSpec(int nLanes, std::vector con // add the channel for the GRP information to the first processor for (auto& spec : pipelines) { o2::tpc::VDriftHelper::requestCCDBInputs(spec.inputs); // add the same CCDB request to each pipeline + if (distortionType) { + spec.inputs.emplace_back("tpcdistortions", o2::header::gDataOriginTPC, "TPCDIST", 0, Lifetime::Condition, ccdbParamSpec(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::DistortionMapMC), {}, 1)); // time-dependent + // load derivative map in case scaling was requested + if (distortionType == 2) { + spec.inputs.emplace_back("tpcdistortionsderiv", o2::header::gDataOriginTPC, "TPCDISTDERIV", 0, Lifetime::Condition, ccdbParamSpec(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::DistortionMapDerivMC), {}, 1)); // time-dependent + } + } } pipelines[0].outputs.emplace_back("TPC", "ROMode", 0, Lifetime::Timeframe); return pipelines; diff --git a/Steer/DigitizerWorkflow/src/TPCDigitizerSpec.h b/Steer/DigitizerWorkflow/src/TPCDigitizerSpec.h index 08a28b2f7647a..65837f7974ba1 100644 --- a/Steer/DigitizerWorkflow/src/TPCDigitizerSpec.h +++ b/Steer/DigitizerWorkflow/src/TPCDigitizerSpec.h @@ -20,9 +20,9 @@ namespace o2 namespace tpc { -o2::framework::DataProcessorSpec getTPCDigitizerSpec(int channel, bool writeGRP, bool mctruth, bool internalwriter); +o2::framework::DataProcessorSpec getTPCDigitizerSpec(int channel, bool writeGRP, bool mctruth, bool internalwriter, int distortionType); -o2::framework::WorkflowSpec getTPCDigitizerSpec(int nLanes, std::vector const& sectors, bool mctruth, bool internalwriter); +o2::framework::WorkflowSpec getTPCDigitizerSpec(int nLanes, std::vector const& sectors, bool mctruth, bool internalwriter, int distortionType); } // end namespace tpc } // end namespace o2 diff --git a/Utilities/DataCompression/test/DataGenerator.h b/Utilities/DataCompression/test/DataGenerator.h index 6c209f9007f6f..3fe0740ccbad8 100644 --- a/Utilities/DataCompression/test/DataGenerator.h +++ b/Utilities/DataCompression/test/DataGenerator.h @@ -154,8 +154,8 @@ class DataGenerator // pointer operator->() const {return &mValue;} // reference operator[](size_type n) const; - bool operator==(const self_type& other) { return mCount == other.mCount; } - bool operator!=(const self_type& other) { return not(*this == other); } + bool operator==(const self_type& other) const { return mCount == other.mCount; } + bool operator!=(const self_type& other) const { return not(*this == other); } private: const ContainerT& mParent; diff --git a/Utilities/ShmManager/src/ShmManager.cxx b/Utilities/ShmManager/src/ShmManager.cxx index 1c54ffcb1552a..866517248277e 100644 --- a/Utilities/ShmManager/src/ShmManager.cxx +++ b/Utilities/ShmManager/src/ShmManager.cxx @@ -152,7 +152,7 @@ struct ShmManager { } #endif - auto ret = regions.emplace(id, make_unique(shmId, id, size)); + auto ret = regions.emplace(id, make_unique(shmId, cfg)); fair::mq::shmem::UnmanagedRegion& region = *(ret.first->second); LOG(info) << "Created unamanged region " << id << " of size " << region.GetSize() << ", starting at " << region.GetData() << ". Locking..."; region.Lock(); diff --git a/prodtests/full-system-test/aggregator-workflow.sh b/prodtests/full-system-test/aggregator-workflow.sh index 9542633bd2cb8..0fc006d183a4d 100755 --- a/prodtests/full-system-test/aggregator-workflow.sh +++ b/prodtests/full-system-test/aggregator-workflow.sh @@ -315,10 +315,6 @@ fi # Forward detectors if [[ $AGGREGATOR_TASKS == FORWARD_TF || $AGGREGATOR_TASKS == ALL ]]; then - # ZDC - if [[ $CALIB_ZDC_TDC == 1 ]]; then - add_W o2-zdc-tdccalib-workflow "" "CalibParamZDC.outputDir=$CALIB_DIR;CalibParamZDC.metaFileDir=$EPN2EOS_METAFILES_DIR" - fi # FT0 if [[ $CALIB_FT0_TIMEOFFSET == 1 ]]; then add_W o2-calibration-ft0-time-offset-calib "--tf-per-slot $FT0_TIMEOFFSET_TF_PER_SLOT --max-delay 0" "FT0CalibParam.mNExtraSlots=0;FT0CalibParam.mRebinFactorPerChID[180]=4;" @@ -338,6 +334,10 @@ if [[ $AGGREGATOR_TASKS == FORWARD_SPORADIC || $AGGREGATOR_TASKS == ALL ]]; then if [[ $CALIB_FDD_INTEGRATEDCURR == 1 ]]; then add_W o2-fdd-merge-integrate-cluster-workflow "--tf-per-slot $INTEGRATEDCURR_TF_PER_SLOT" fi + # ZDC + if [[ $CALIB_ZDC_TDC == 1 ]]; then + add_W o2-zdc-tdccalib-workflow "" "CalibParamZDC.outputDir=$CALIB_DIR;CalibParamZDC.metaFileDir=$EPN2EOS_METAFILES_DIR" + fi fi if [[ "${GEN_TOPO_VERBOSE:-}" == "1" ]]; then diff --git a/prodtests/full-system-test/dpl-workflow.sh b/prodtests/full-system-test/dpl-workflow.sh index a5416204e27fa..2a03a0ebc12f0 100755 --- a/prodtests/full-system-test/dpl-workflow.sh +++ b/prodtests/full-system-test/dpl-workflow.sh @@ -272,10 +272,9 @@ GPU_CONFIG_SELF="--severity $SEVERITY_TPC" parse_TPC_CORR_SCALING() { while [[ $# -gt 0 ]]; do - echo $1 case "$1" in - --lumi-type=*) TPC_CORR_OPT+=" --lumi-type ${1#*=}"; shift 1;; - --lumi-type) TPC_CORR_OPT+=" --lumi-type ${2}"; shift 2;; + --lumi-type=*) TPC_CORR_OPT+=" --lumi-type ${1#*=}"; [[ ${1#*=} == "2" ]] && NEED_TPC_SCALERS_WF=1; shift 1;; + --lumi-type) TPC_CORR_OPT+=" --lumi-type ${2}"; [[ ${2} == "2" ]] && NEED_TPC_SCALERS_WF=1; shift 2;; --corrmap-lumi-mode=*) TPC_CORR_OPT+=" --corrmap-lumi-mode ${1#*=}"; shift 1;; --corrmap-lumi-mode) TPC_CORR_OPT+=" --corrmap-lumi-mode ${2}"; shift 2;; *) TPC_CORR_KEY+="$1;"; shift 1;; diff --git a/prodtests/full-system-test/shm-tool.sh b/prodtests/full-system-test/shm-tool.sh index 9dc5c7de7137b..8c7ac08f82dbc 100755 --- a/prodtests/full-system-test/shm-tool.sh +++ b/prodtests/full-system-test/shm-tool.sh @@ -1,4 +1,8 @@ #!/bin/bash +if [[ -z $SHM_MANAGER_SHMID ]]; then + echo You must provide SHM_MANAGER_SHMID + exit 1 +fi if [[ "0$EPN_NODE_MI100" == "01" && -z $EPN_GLOBAL_SCALING ]]; then EPN_GLOBAL_SCALING="3 / 2" diff --git a/prodtests/full-system-test/start_tmux.sh b/prodtests/full-system-test/start_tmux.sh index f45a8a83841ec..9c1a6e1a473ef 100755 --- a/prodtests/full-system-test/start_tmux.sh +++ b/prodtests/full-system-test/start_tmux.sh @@ -78,10 +78,10 @@ if [ $1 == "dd" ]; then export GPU_NUM_MEM_REG_CALLBACKS=$(($NUM_DPL_WORKFLOWS + 3)) elif [ $1 == "tf" ]; then export CMD=tf-reader.sh - export GPU_NUM_MEM_REG_CALLBACKS=$((NUM_DPL_WORKFLOWS + 1)) + export GPU_NUM_MEM_REG_CALLBACKS=$((NUM_DPL_WORKFLOWS + ${NUMAGPUIDS:-0})) elif [ $1 == "rr" ]; then export CMD=raw-reader.sh - export GPU_NUM_MEM_REG_CALLBACKS=$(($NUM_DPL_WORKFLOWS + 1)) + export GPU_NUM_MEM_REG_CALLBACKS=$(($NUM_DPL_WORKFLOWS + ${NUMAGPUIDS:-0})) fi if [ "0$FST_TMUX_NOWAIT" != "01" ]; then @@ -112,11 +112,6 @@ else fi [[ ! -z $FST_TMUX_DD_WAIT ]] && FST_SLEEP2=$FST_TMUX_DD_WAIT -if [[ ! -z $FST_TMUX_SINGLENUMA ]]; then - eval "FST_SLEEP$((FST_TMUX_SINGLENUMA ^ 1))=\"0; echo SKIPPED; sleep 1000; exit\"" - export GPU_NUM_MEM_REG_CALLBACKS=$(($GPU_NUM_MEM_REG_CALLBACKS - 1)) -fi - if workflow_has_parameter CALIB_PROXIES; then CALIB_COMMAND="$GEN_TOPO_MYDIR/aggregator-workflow.sh" if [[ -z ${CALIB_TASKS:-} ]]; then diff --git a/run/O2SimDevice.h b/run/O2SimDevice.h index ef1812b11d278..09ff0f5fe97d3 100644 --- a/run/O2SimDevice.h +++ b/run/O2SimDevice.h @@ -23,7 +23,6 @@ #include "TVirtualMC.h" #include "TMessage.h" #include -#include #include #include #include @@ -258,7 +257,6 @@ class O2SimDevice final : public fair::mq::Device << "part " << info.part << "/" << info.nparts; LOG(info) << workerStr() << " Setting seed for this sub-event to " << chunk->mSubEventInfo.seed; gRandom->SetSeed(chunk->mSubEventInfo.seed); - o2::base::VMCSeederService::instance().setSeed(); // Process one event auto& conf = o2::conf::SimConfig::Instance();