diff --git a/.github/workflows/mega-linter.yml b/.github/workflows/mega-linter.yml
index 0b700b2a01d..e08d00c2cc2 100644
--- a/.github/workflows/mega-linter.yml
+++ b/.github/workflows/mega-linter.yml
@@ -1,6 +1,6 @@
---
# MegaLinter GitHub Action configuration file
-# More info at https://oxsecurity.github.io/megalinter
+# More info at https://megalinter.io
name: MegaLinter
'on': [pull_request_target]
@@ -17,7 +17,7 @@ concurrency:
cancel-in-progress: true
jobs:
- build:
+ megalinter:
name: MegaLinter
runs-on: ubuntu-latest
steps:
@@ -38,11 +38,11 @@ jobs:
- name: MegaLinter
id: ml
# You can override MegaLinter flavor used to have faster performances
- # More info at https://oxsecurity.github.io/megalinter/flavors/
- uses: oxsecurity/megalinter@v6
+ # More info at https://megalinter.io/flavors/
+ uses: oxsecurity/megalinter@v7
env:
# All available variables are described in documentation:
- # https://oxsecurity.github.io/megalinter/configuration/
+ # https://megalinter.io/configuration/
# Validates all source when push on master, else just the diff with
# master. Override with true if you always want to lint all sources.
VALIDATE_ALL_CODEBASE: false
@@ -83,7 +83,7 @@ jobs:
yourself and update the pull request, or merge this PR in yours.
You can find how to run MegaLinter locally at
- .
+ .
# We do not create PRs if the branch is not there.
continue-on-error: true
diff --git a/.gitignore b/.gitignore
index 14be3e78db7..4bc463b2cfc 100644
--- a/.gitignore
+++ b/.gitignore
@@ -7,3 +7,4 @@ compile_commands.json
megalinter-reports/
.mypy_cache/
github_conf/
+.cache
diff --git a/.mega-linter.yml b/.mega-linter.yml
index 9e3581c67c9..7706edad085 100644
--- a/.mega-linter.yml
+++ b/.mega-linter.yml
@@ -17,10 +17,11 @@ DISABLE_LINTERS:
- YAML_V8R
- YAML_PRETTIER
- REPOSITORY_DEVSKIM
+ - REPOSITORY_KICS
- REPOSITORY_SECRETLINT
- REPOSITORY_TRIVY
DISABLE_ERRORS_LINTERS: # If errors are found by these linters, they will be considered as non blocking.
- - PYTHON_BANDIT # The bandit check is overly broad and complains about subprocess usage in Scripts/update_ccdb.py.
+ - PYTHON_BANDIT # The bandit check is overly broad and complains about subprocess usage.
SHOW_ELAPSED_TIME: true
FILEIO_REPORTER: false
GITHUB_COMMENT_REPORTER: false
@@ -28,3 +29,6 @@ UPDATED_SOURCES_REPORTER: false
PRINT_ALPACA: false # Don't print ASCII alpaca in the log
PRINT_ALL_FILES: true # Print all processed files
FLAVOR_SUGGESTIONS: false # Don't show suggestions about different MegaLinter flavors
+PYTHON_ISORT_CONFIG_FILE: pyproject.toml
+PYTHON_PYRIGHT_CONFIG_FILE: pyproject.toml
+PYTHON_RUFF_CONFIG_FILE: pyproject.toml
diff --git a/ALICE3/Core/DelphesO2TrackSmearer.cxx b/ALICE3/Core/DelphesO2TrackSmearer.cxx
index 74e14591f74..f4ea4ebdf3c 100644
--- a/ALICE3/Core/DelphesO2TrackSmearer.cxx
+++ b/ALICE3/Core/DelphesO2TrackSmearer.cxx
@@ -110,7 +110,7 @@ bool TrackSmearer::loadTable(int pdg, const char* filename, bool forceReload)
/*****************************************************************/
lutEntry_t*
- TrackSmearer::getLUTEntry(int pdg, float nch, float radius, float eta, float pt)
+ TrackSmearer::getLUTEntry(int pdg, float nch, float radius, float eta, float pt, float& interpolatedEff)
{
auto ipdg = getIndexPDG(pdg);
if (!mLUTHeader[ipdg])
@@ -119,13 +119,56 @@ lutEntry_t*
auto irad = mLUTHeader[ipdg]->radmap.find(radius);
auto ieta = mLUTHeader[ipdg]->etamap.find(eta);
auto ipt = mLUTHeader[ipdg]->ptmap.find(pt);
+
+ // Interpolate if requested
+ auto fraction = mLUTHeader[ipdg]->nchmap.fracPositionWithinBin(nch);
+ if (mInterpolateEfficiency) {
+ if (fraction > 0.5) {
+ if (mWhatEfficiency == 1) {
+ if (inch < mLUTHeader[ipdg]->nchmap.nbins - 1) {
+ interpolatedEff = (1.5f - fraction) * mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff + (-0.5f + fraction) * mLUTEntry[ipdg][inch + 1][irad][ieta][ipt]->eff;
+ } else {
+ interpolatedEff = mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff;
+ }
+ }
+ if (mWhatEfficiency == 2) {
+ if (inch < mLUTHeader[ipdg]->nchmap.nbins - 1) {
+ interpolatedEff = (1.5f - fraction) * mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff2 + (-0.5f + fraction) * mLUTEntry[ipdg][inch + 1][irad][ieta][ipt]->eff2;
+ } else {
+ interpolatedEff = mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff2;
+ }
+ }
+ } else {
+ float comparisonValue = mLUTHeader[ipdg]->nchmap.log ? log10(nch) : nch;
+ if (mWhatEfficiency == 1) {
+ if (inch > 0 && comparisonValue < mLUTHeader[ipdg]->nchmap.max) {
+ interpolatedEff = (0.5f + fraction) * mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff + (0.5f - fraction) * mLUTEntry[ipdg][inch - 1][irad][ieta][ipt]->eff;
+ } else {
+ interpolatedEff = mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff;
+ }
+ }
+ if (mWhatEfficiency == 2) {
+ if (inch > 0 && comparisonValue < mLUTHeader[ipdg]->nchmap.max) {
+ interpolatedEff = (0.5f + fraction) * mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff2 + (0.5f - fraction) * mLUTEntry[ipdg][inch - 1][irad][ieta][ipt]->eff2;
+ } else {
+ interpolatedEff = mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff2;
+ }
+ }
+ }
+ } else {
+ if (mWhatEfficiency == 1)
+ interpolatedEff = mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff;
+ if (mWhatEfficiency == 2)
+ interpolatedEff = mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff2;
+ }
return mLUTEntry[ipdg][inch][irad][ieta][ipt];
} //;
/*****************************************************************/
-bool TrackSmearer::smearTrack(O2Track& o2track, lutEntry_t* lutEntry)
+bool TrackSmearer::smearTrack(O2Track& o2track, lutEntry_t* lutEntry, float interpolatedEff)
{
+ bool isReconstructed = true;
// generate efficiency
if (mUseEfficiency) {
auto eff = 0.;
@@ -133,9 +176,16 @@ bool TrackSmearer::smearTrack(O2Track& o2track, lutEntry_t* lutEntry)
eff = lutEntry->eff;
if (mWhatEfficiency == 2)
eff = lutEntry->eff2;
+ if (mInterpolateEfficiency)
+ eff = interpolatedEff;
if (gRandom->Uniform() > eff)
- return false;
+ isReconstructed = false;
}
+
+ // return false already now in case not reco'ed
+ if (!isReconstructed && mSkipUnreconstructed)
+ return false;
+
// transform params vector and smear
double params_[5];
for (int i = 0; i < 5; ++i) {
@@ -158,7 +208,7 @@ bool TrackSmearer::smearTrack(O2Track& o2track, lutEntry_t* lutEntry)
// set covariance matrix
for (int i = 0; i < 15; ++i)
o2track.setCov(lutEntry->covm[i], i);
- return true;
+ return isReconstructed;
}
/*****************************************************************/
@@ -171,31 +221,62 @@ bool TrackSmearer::smearTrack(O2Track& o2track, int pdg, float nch)
pt *= 2.f;
}
auto eta = o2track.getEta();
- auto lutEntry = getLUTEntry(pdg, nch, 0., eta, pt);
+ float interpolatedEff = 0.0f;
+ auto lutEntry = getLUTEntry(pdg, nch, 0., eta, pt, interpolatedEff);
if (!lutEntry || !lutEntry->valid)
return false;
- return smearTrack(o2track, lutEntry);
+ return smearTrack(o2track, lutEntry, interpolatedEff);
}
/*****************************************************************/
// relative uncertainty on pt
double TrackSmearer::getPtRes(int pdg, float nch, float eta, float pt)
{
- auto lutEntry = getLUTEntry(pdg, nch, 0., eta, pt);
+ float dummy = 0.0f;
+ auto lutEntry = getLUTEntry(pdg, nch, 0., eta, pt, dummy);
auto val = sqrt(lutEntry->covm[14]) * lutEntry->pt;
return val;
}
+
/*****************************************************************/
// relative uncertainty on eta
double TrackSmearer::getEtaRes(int pdg, float nch, float eta, float pt)
{
- auto lutEntry = getLUTEntry(pdg, nch, 0., eta, pt);
- auto sigmatgl = sqrt(lutEntry->covm[9]); // sigmatgl2
- auto etaRes = 1 / (sqrt(1 + sigmatgl * sigmatgl)); // propagate tgl to eta uncertainty
- etaRes /= lutEntry->eta; // relative uncertainty
+ float dummy = 0.0f;
+ auto lutEntry = getLUTEntry(pdg, nch, 0., eta, pt, dummy);
+ auto sigmatgl = sqrt(lutEntry->covm[9]); // sigmatgl2
+ auto etaRes = fabs(sin(2.0 * atan(exp(-eta)))) * sigmatgl; // propagate tgl to eta uncertainty
+ etaRes /= lutEntry->eta; // relative uncertainty
return etaRes;
}
+/*****************************************************************/
+// absolute uncertainty on pt
+double TrackSmearer::getAbsPtRes(int pdg, float nch, float eta, float pt)
+{
+ float dummy = 0.0f;
+ auto lutEntry = getLUTEntry(pdg, nch, 0., eta, pt, dummy);
+ auto val = sqrt(lutEntry->covm[14]) * pow(lutEntry->pt, 2);
+ return val;
+}
+/*****************************************************************/
+// absolute uncertainty on eta
+double TrackSmearer::getAbsEtaRes(int pdg, float nch, float eta, float pt)
+{
+ float dummy = 0.0f;
+ auto lutEntry = getLUTEntry(pdg, nch, 0., eta, pt, dummy);
+ auto sigmatgl = sqrt(lutEntry->covm[9]); // sigmatgl2
+ auto etaRes = fabs(sin(2.0 * atan(exp(-eta)))) * sigmatgl; // propagate tgl to eta uncertainty
+ return etaRes;
+}
+/*****************************************************************/
+// efficiency
+double TrackSmearer::getEfficiency(int pdg, float nch, float eta, float pt)
+{
+ float efficiency = 0.0f;
+ getLUTEntry(pdg, nch, 0., eta, pt, efficiency);
+ return efficiency;
+}
/*****************************************************************/
// Only in DelphesO2
// bool TrackSmearer::smearTrack(Track& track, bool atDCA)
diff --git a/ALICE3/Core/DelphesO2TrackSmearer.h b/ALICE3/Core/DelphesO2TrackSmearer.h
index 279fd9e3046..22a0bbd1346 100644
--- a/ALICE3/Core/DelphesO2TrackSmearer.h
+++ b/ALICE3/Core/DelphesO2TrackSmearer.h
@@ -53,7 +53,23 @@ struct map_t {
if (log)
return pow(10., val);
return val;
- } //;
+ }
+ // function needed to interpolate some dimensions
+ float fracPositionWithinBin(float val)
+ {
+ float width = (max - min) / nbins;
+ int bin;
+ float returnVal = 0.5f;
+ if (log) {
+ bin = static_cast((log10(val) - min) / width);
+ returnVal = ((log10(val) - min) / width) - bin;
+ } else {
+ bin = static_cast((val - min) / width);
+ returnVal = val / width - bin;
+ }
+ return returnVal;
+ }
+
int find(float val)
{
float width = (max - min) / nbins;
@@ -162,15 +178,20 @@ class TrackSmearer
/** LUT methods **/
bool loadTable(int pdg, const char* filename, bool forceReload = false);
void useEfficiency(bool val) { mUseEfficiency = val; } //;
+ void interpolateEfficiency(bool val) { mInterpolateEfficiency = val; } //;
+ void skipUnreconstructed(bool val) { mSkipUnreconstructed = val; } //;
void setWhatEfficiency(int val) { mWhatEfficiency = val; } //;
lutHeader_t* getLUTHeader(int pdg) { return mLUTHeader[getIndexPDG(pdg)]; } //;
- lutEntry_t* getLUTEntry(int pdg, float nch, float radius, float eta, float pt);
+ lutEntry_t* getLUTEntry(int pdg, float nch, float radius, float eta, float pt, float& interpolatedEff);
- bool smearTrack(O2Track& o2track, lutEntry_t* lutEntry);
+ bool smearTrack(O2Track& o2track, lutEntry_t* lutEntry, float interpolatedEff);
bool smearTrack(O2Track& o2track, int pdg, float nch);
// bool smearTrack(Track& track, bool atDCA = true); // Only in DelphesO2
double getPtRes(int pdg, float nch, float eta, float pt);
double getEtaRes(int pdg, float nch, float eta, float pt);
+ double getAbsPtRes(int pdg, float nch, float eta, float pt);
+ double getAbsEtaRes(int pdg, float nch, float eta, float pt);
+ double getEfficiency(int pdg, float nch, float eta, float pt);
int getIndexPDG(int pdg)
{
@@ -203,6 +224,8 @@ class TrackSmearer
lutHeader_t* mLUTHeader[nLUTs] = {nullptr};
lutEntry_t***** mLUTEntry[nLUTs] = {nullptr};
bool mUseEfficiency = true;
+ bool mInterpolateEfficiency = false;
+ bool mSkipUnreconstructed = true; // don't smear tracks that are not reco'ed
int mWhatEfficiency = 1;
float mdNdEta = 1600.;
};
diff --git a/ALICE3/DataModel/A3DecayFinderTables.h b/ALICE3/DataModel/A3DecayFinderTables.h
new file mode 100644
index 00000000000..cc7e8f3d3c7
--- /dev/null
+++ b/ALICE3/DataModel/A3DecayFinderTables.h
@@ -0,0 +1,61 @@
+// 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 A3DecayFinderTables.h
+/// \since 04/07/2023
+/// \brief Set of tables for ALICE 3 decay finder
+///
+
+#ifndef ALICE3_DATAMODEL_A3DECAYFINDERTABLES_H_
+#define ALICE3_DATAMODEL_A3DECAYFINDERTABLES_H_
+
+// O2 includes
+#include "Framework/AnalysisDataModel.h"
+
+enum a3selectionBit : uint32_t { kDCAxy = 0,
+ kInnerTOFPion,
+ kInnerTOFKaon,
+ kInnerTOFProton,
+ kOuterTOFPion,
+ kOuterTOFKaon,
+ kOuterTOFProton,
+ kRICHPion,
+ kRICHKaon,
+ kRICHProton,
+ kTruePion,
+ kTrueKaon,
+ kTrueProton,
+ kTruePiPlusFromD,
+ kTrueKaPlusFromD,
+ kTruePiMinusFromD,
+ kTrueKaMinusFromD,
+ kTruePiPlusFromLc,
+ kTrueKaPlusFromLc,
+ kTruePrPlusFromLc,
+ kTruePiMinusFromLc,
+ kTrueKaMinusFromLc,
+ kTruePrMinusFromLc };
+
+namespace o2::aod
+{
+namespace a3DecayMap
+{
+DECLARE_SOA_COLUMN(DecayMap, decayMap, uint32_t); //! simple map to process passing / not passing criteria
+} // namespace a3DecayMap
+DECLARE_SOA_TABLE(Alice3DecayMaps, "AOD", "ALICE3DECAYMAP",
+ a3DecayMap::DecayMap);
+
+using Alice3DecayMap = Alice3DecayMaps::iterator;
+
+} // namespace o2::aod
+
+#endif // ALICE3_DATAMODEL_A3DECAYFINDERTABLES_H_
diff --git a/ALICE3/DataModel/OTFRICH.h b/ALICE3/DataModel/OTFRICH.h
new file mode 100644
index 00000000000..e853a2ff30f
--- /dev/null
+++ b/ALICE3/DataModel/OTFRICH.h
@@ -0,0 +1,46 @@
+// 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 OTFRICH.h
+/// \author Nicola Nicassio, University and INFN Bari
+/// \since 15/05/2023
+/// \brief Set of tables for the ALICE 3 OTFRICH information
+///
+
+#ifndef ALICE3_DATAMODEL_OTFRICH_H_
+#define ALICE3_DATAMODEL_OTFRICH_H_
+
+// O2 includes
+#include "Framework/AnalysisDataModel.h"
+
+namespace o2::aod
+{
+namespace upgrade_rich
+{
+DECLARE_SOA_COLUMN(NSigmaElectronRich, nSigmaElectronRich, float); //! NSigma electron BarrelRich
+DECLARE_SOA_COLUMN(NSigmaMuonRich, nSigmaMuonRich, float); //! NSigma muon BarrelRich
+DECLARE_SOA_COLUMN(NSigmaPionRich, nSigmaPionRich, float); //! NSigma pion BarrelRich
+DECLARE_SOA_COLUMN(NSigmaKaonRich, nSigmaKaonRich, float); //! NSigma kaon BarrelRich
+DECLARE_SOA_COLUMN(NSigmaProtonRich, nSigmaProtonRich, float); //! NSigma proton BarrelRich
+} // namespace upgrade_rich
+DECLARE_SOA_TABLE(UpgradeRichs, "AOD", "UPGRADERICH",
+ upgrade_rich::NSigmaElectronRich,
+ upgrade_rich::NSigmaMuonRich,
+ upgrade_rich::NSigmaPionRich,
+ upgrade_rich::NSigmaKaonRich,
+ upgrade_rich::NSigmaProtonRich);
+
+using UpgradeRich = UpgradeRichs::iterator;
+
+} // namespace o2::aod
+
+#endif // ALICE3_DATAMODEL_OTFRICH_H_
diff --git a/ALICE3/DataModel/tracksAlice3.h b/ALICE3/DataModel/tracksAlice3.h
new file mode 100644
index 00000000000..c3a3e82a994
--- /dev/null
+++ b/ALICE3/DataModel/tracksAlice3.h
@@ -0,0 +1,38 @@
+// 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 tracksAlice3.h
+/// \author David Dobrigkeit Chinellato
+/// \since 11/05/2023
+/// \brief Table for ALICE 3 track-related info
+///
+
+#ifndef ALICE3_DATAMODEL_TRACKSALICE3_H_
+#define ALICE3_DATAMODEL_TRACKSALICE3_H_
+
+// O2 includes
+#include "Framework/AnalysisDataModel.h"
+
+namespace o2::aod
+{
+namespace track_alice3
+{
+DECLARE_SOA_COLUMN(IsReconstructed, isReconstructed, bool); //! is reconstructed or not
+} // namespace track_alice3
+DECLARE_SOA_TABLE(TracksAlice3, "AOD", "TRACKSALICE3",
+ track_alice3::IsReconstructed);
+
+using TrackAlice3 = TracksAlice3::iterator;
+
+} // namespace o2::aod
+
+#endif // ALICE3_DATAMODEL_TRACKSALICE3_H_
diff --git a/ALICE3/TableProducer/CMakeLists.txt b/ALICE3/TableProducer/CMakeLists.txt
index 72178be3179..52300bfc644 100644
--- a/ALICE3/TableProducer/CMakeLists.txt
+++ b/ALICE3/TableProducer/CMakeLists.txt
@@ -9,6 +9,8 @@
# granted to it by virtue of its status as an Intergovernmental Organization
# or submit itself to any jurisdiction.
+add_subdirectory(OTF)
+
o2physics_add_dpl_workflow(alice3-trackselection
SOURCES alice3-trackselection.cxx
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore
@@ -29,12 +31,8 @@ o2physics_add_dpl_workflow(alice3-centrality
PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore
COMPONENT_NAME Analysis)
-o2physics_add_dpl_workflow(onthefly-tracker
- SOURCES onTheFlyTracker.cxx
- PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2::ReconstructionDataFormats O2::DetectorsCommonDataFormats O2Physics::ALICE3Core
+o2physics_add_dpl_workflow(alice3-decayfinder
+ SOURCES alice3-decayfinder.cxx
+ PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DCAFitter
COMPONENT_NAME Analysis)
-o2physics_add_dpl_workflow(onthefly-tofpid
- SOURCES onTheFlyTOFPID.cxx
- PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2::ReconstructionDataFormats O2::DetectorsCommonDataFormats
- COMPONENT_NAME Analysis)
diff --git a/ALICE3/TableProducer/OTF/CMakeLists.txt b/ALICE3/TableProducer/OTF/CMakeLists.txt
new file mode 100644
index 00000000000..305fcce4b66
--- /dev/null
+++ b/ALICE3/TableProducer/OTF/CMakeLists.txt
@@ -0,0 +1,25 @@
+# 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.
+
+o2physics_add_dpl_workflow(onthefly-tracker
+ SOURCES onTheFlyTracker.cxx
+ PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2::ReconstructionDataFormats O2::DetectorsCommonDataFormats O2::DetectorsVertexing O2Physics::ALICE3Core
+ COMPONENT_NAME Analysis)
+
+o2physics_add_dpl_workflow(onthefly-tofpid
+ SOURCES onTheFlyTOFPID.cxx
+ PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2::ReconstructionDataFormats O2::DetectorsCommonDataFormats O2Physics::ALICE3Core
+ COMPONENT_NAME Analysis)
+
+o2physics_add_dpl_workflow(onthefly-richpid
+ SOURCES onTheFlyRICHPID.cxx
+ PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2::ReconstructionDataFormats O2::DetectorsCommonDataFormats O2Physics::ALICE3Core
+ COMPONENT_NAME Analysis)
diff --git a/ALICE3/TableProducer/OTF/onTheFlyRICHPID.cxx b/ALICE3/TableProducer/OTF/onTheFlyRICHPID.cxx
new file mode 100644
index 00000000000..82090899472
--- /dev/null
+++ b/ALICE3/TableProducer/OTF/onTheFlyRICHPID.cxx
@@ -0,0 +1,624 @@
+// 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.
+
+//
+// Task to add a table of track parameters propagated to the primary vertex
+//
+
+#include
+#include
+#include "Framework/AnalysisDataModel.h"
+#include "Framework/AnalysisTask.h"
+#include "Framework/runDataProcessing.h"
+#include "Framework/RunningWorkflowInfo.h"
+#include "Framework/HistogramRegistry.h"
+#include "Framework/O2DatabasePDGPlugin.h"
+#include "Framework/ASoAHelpers.h"
+#include "Common/DataModel/TrackSelectionTables.h"
+#include "Common/Core/trackUtilities.h"
+#include "ReconstructionDataFormats/DCA.h"
+#include "DetectorsBase/Propagator.h"
+#include "DetectorsBase/GeometryManager.h"
+#include "CommonUtils/NameConf.h"
+#include "CCDB/CcdbApi.h"
+#include "CCDB/BasicCCDBManager.h"
+#include "DataFormatsParameters/GRPMagField.h"
+#include "DataFormatsCalibration/MeanVertexObject.h"
+#include "CommonConstants/GeomConstants.h"
+#include "CommonConstants/PhysicsConstants.h"
+#include "TRandom3.h"
+#include "TVector3.h"
+#include "TString.h"
+#include "ALICE3/DataModel/OTFRICH.h"
+#include "DetectorsVertexing/HelixHelper.h"
+
+#include "TableHelper.h"
+#include "ALICE3/Core/DelphesO2TrackSmearer.h"
+
+/// \file onTheFlyRichPid.cxx
+///
+/// \brief This task goes straight from a combination of track table and mcParticles
+/// and a projective bRICH configuration to a table of TOF NSigmas for the particles
+/// being analysed. It currently contemplates 5 particle types:
+/// electrons, pions, kaons, protons and muons.
+///
+/// More particles could be added but would have to be added to the LUT
+/// being used in the onTheFly tracker task.
+///
+/// \warning Geometry parameters are configurable, but resolution values should be adapted.
+/// Since angular resolution depends on the specific geometric details, it is better to
+/// calculate it from full simulation and add new input. Alternatively, an analytical
+/// expression can be provided as a function of the main parameters.
+///
+/// \author David Dobrigkeit Chinellato, UNICAMP, Nicola Nicassio, University and INFN Bari
+
+using namespace o2;
+using namespace o2::framework;
+
+struct OnTheFlyRichPid {
+ Produces upgradeRich;
+
+ // necessary for particle charges
+ Service pdg;
+
+ // master setting: magnetic field
+ Configurable dBz{"dBz", 20, "magnetic field (kilogauss)"};
+
+ // add rich-specific configurables here
+ Configurable bRichNumberOfSectors{"bRichNumberOfSectors", 21, "barrel RICH number of sectors"};
+ Configurable bRichPhotodetectorCentralModuleHalfLength{"bRichPhotodetectorCentralModuleHalfLength", 18.4 / 2.0, "barrel RICH photodetector central module half length (cm)"};
+ Configurable bRichPhotodetectorOtherModuleLength{"bRichPhotodetectorOtherModuleLength", 18.4, "barrel RICH photodetector other module length (cm)"};
+ Configurable bRichRadiatorInnerRadius{"bRichRadiatorInnerRadius", 85., "barrel RICH radiator inner radius (cm)"};
+ Configurable bRichPhotodetectorOuterRadius{"bRichPhotodetectorOuterRadius", 112., "barrel RICH photodetector outer radius (cm)"};
+ Configurable bRichRadiatorThickness{"bRichRadiatorThickness", 2., "barrel RICH radiator thickness (cm)"};
+ Configurable bRichRefractiveIndex{"bRichRefractiveIndex", 1.03, "barrel RICH refractive index"};
+ Configurable bRichFlagAbsorbingWalls{"bRichFlagAbsorbingWalls", false, "barrel RICH flag absorbing walls between sectors"};
+ Configurable nStepsLIntegrator{"nStepsLIntegrator", 200, "number of steps in length integrator"};
+ Configurable doQAplots{"doQAplots", true, "do basic velocity plot qa"};
+ Configurable nBinsThetaRing{"nBinsThetaRing", 3000, "number of bins in theta ring"};
+ Configurable nBinsP{"nBinsP", 400, "number of bins in momentum"};
+ Configurable nBinsNsigmaCorrectSpecies{"nBinsNsigmaCorrectSpecies", 200, "number of bins in Nsigma plot (correct speies)"};
+ Configurable nBinsNsigmaWrongSpecies{"nBinsNsigmaWrongSpecies", 200, "number of bins in Nsigma plot (wrong species)"};
+ Configurable nBinsAngularRes{"nBinsAngularRes", 400, "number of bins plots angular resolution"};
+ Configurable nBinsRelativeEtaPt{"nBinsRelativeEtaPt", 400, "number of bins plots pt and eta relative errors"};
+ Configurable nBinsEta{"nBinsEta", 400, "number of bins plot relative eta error"};
+ Configurable flagIncludeTrackAngularRes{"flagIncludeTrackAngularRes", true, "flag to include or exclude track time resolution"};
+ Configurable multiplicityEtaRange{"multiplicityEtaRange", 0.800000012, "eta range to compute the multiplicity"};
+ Configurable flagRICHLoadDelphesLUTs{"flagRICHLoadDelphesLUTs", false, "flag to load Delphes LUTs for tracking correction (use recoTrack parameters if false)"};
+
+ Configurable lutEl{"lutEl", "lutCovm.el.dat", "LUT for electrons"};
+ Configurable lutMu{"lutMu", "lutCovm.mu.dat", "LUT for muons"};
+ Configurable lutPi{"lutPi", "lutCovm.pi.dat", "LUT for pions"};
+ Configurable lutKa{"lutKa", "lutCovm.ka.dat", "LUT for kaons"};
+ Configurable lutPr{"lutPr", "lutCovm.pr.dat", "LUT for protons"};
+
+ o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrNONE;
+
+ // Track smearer (here used to get relative pt and eta uncertainties)
+ o2::delphes::DelphesO2TrackSmearer mSmearer;
+
+ // needed: random number generator for smearing
+ TRandom3 pRandomNumberGenerator;
+
+ // for handling basic QA histograms if requested
+ HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject};
+
+ /// Flag unphysical and unavailable values
+ float error_value = -1000;
+
+ // Variables projective/hybrid layout
+ int mNumberSectors = bRichNumberOfSectors; // 21;
+ float mTileLength = bRichPhotodetectorOtherModuleLength; // 18.4; // [cm]
+ float mTileLengthCentral = bRichPhotodetectorCentralModuleHalfLength; // 18.4 / 2.0; // [cm]
+ float mProjectiveLengthInner = mTileLengthCentral;
+ float mRadiusProjIn = bRichRadiatorInnerRadius; // 85.; // [cm]
+ float mRadiusProjOut = bRichPhotodetectorOuterRadius; // 112.; // [cm]
+ std::vector det_centers;
+ std::vector rad_centers;
+ std::vector angle_centers;
+ std::vector theta_min;
+ std::vector theta_max;
+
+ // Update projective geometry
+ void updateProjectiveParameters()
+ {
+ const int number_of_sectors_in_z = mNumberSectors;
+ det_centers.resize(number_of_sectors_in_z);
+ rad_centers.resize(number_of_sectors_in_z);
+ angle_centers.resize(number_of_sectors_in_z);
+ theta_min.resize(number_of_sectors_in_z);
+ theta_max.resize(number_of_sectors_in_z);
+ float square_size_barrel_cylinder = 2.0 * mTileLengthCentral;
+ float square_size_z = mTileLength;
+ float R_min = mRadiusProjIn;
+ float R_max = mRadiusProjOut;
+ std::vector theta_bi;
+ std::vector R0_tilt;
+ std::vector z0_tilt;
+ std::vector T_r_plus_g;
+ std::vector l_aerogel_z;
+ std::vector l_detector_z;
+ theta_bi.resize(number_of_sectors_in_z);
+ R0_tilt.resize(number_of_sectors_in_z);
+ z0_tilt.resize(number_of_sectors_in_z);
+ T_r_plus_g.resize(number_of_sectors_in_z);
+ l_aerogel_z.resize(number_of_sectors_in_z);
+ l_detector_z.resize(number_of_sectors_in_z);
+
+ // Central sector
+ int i_central_mirror = static_cast((number_of_sectors_in_z) / 2.0);
+ float m_val = std::tan(0.0);
+ theta_bi[i_central_mirror] = std::atan(m_val);
+ R0_tilt[i_central_mirror] = R_max;
+ z0_tilt[i_central_mirror] = R0_tilt[i_central_mirror] * std::tan(theta_bi[i_central_mirror]);
+ l_detector_z[i_central_mirror] = square_size_barrel_cylinder;
+ l_aerogel_z[i_central_mirror] = std::sqrt(1.0 + m_val * m_val) * R_min * square_size_barrel_cylinder / (std::sqrt(1.0 + m_val * m_val) * R_max - m_val * square_size_barrel_cylinder);
+ T_r_plus_g[i_central_mirror] = R_max - R_min;
+ float t = std::tan(atan(m_val) + std::atan(square_size_barrel_cylinder / (2.0 * R_max * std::sqrt(1.0 + m_val * m_val) - square_size_barrel_cylinder * m_val)));
+ theta_max[i_central_mirror] = M_PI / 2.0 - std::atan(t);
+ theta_min[i_central_mirror] = M_PI / 2.0 + std::atan(t);
+ mProjectiveLengthInner = R_min * t;
+ for (int i = static_cast((number_of_sectors_in_z) / 2.0) + 1; i < number_of_sectors_in_z; i++) {
+ float par_a = t;
+ float par_b = 2.0 * R_max / square_size_z;
+ m_val = (std::sqrt(par_a * par_a * par_b * par_b + par_b * par_b - 1.0) + par_a * par_b * par_b) / (par_b * par_b - 1.0);
+ theta_min[i] = M_PI / 2.0 - std::atan(t);
+ theta_max[2 * i_central_mirror - i] = M_PI / 2.0 + std::atan(t);
+ t = std::tan(std::atan(m_val) + std::atan(square_size_z / (2.0 * R_max * std::sqrt(1.0 + m_val * m_val) - square_size_z * m_val)));
+ theta_max[i] = M_PI / 2.0 - std::atan(t);
+ theta_min[2 * i_central_mirror - i] = M_PI / 2.0 + std::atan(t);
+ // Forward sectors
+ theta_bi[i] = std::atan(m_val);
+ R0_tilt[i] = R_max - square_size_z / 2.0 * std::sin(std::atan(m_val));
+ z0_tilt[i] = R0_tilt[i] * std::tan(theta_bi[i]);
+ l_detector_z[i] = square_size_z;
+ l_aerogel_z[i] = std::sqrt(1.0 + m_val * m_val) * R_min * square_size_z / (std::sqrt(1.0 + m_val * m_val) * R_max - m_val * square_size_z);
+ T_r_plus_g[i] = std::sqrt(1.0 + m_val * m_val) * (R_max - R_min) - m_val / 2.0 * (square_size_z + l_aerogel_z[i]);
+ // Backword sectors
+ theta_bi[2 * i_central_mirror - i] = -std::atan(m_val);
+ R0_tilt[2 * i_central_mirror - i] = R_max - square_size_z / 2.0 * std::sin(std::atan(m_val));
+ z0_tilt[2 * i_central_mirror - i] = -R0_tilt[i] * std::tan(theta_bi[i]);
+ l_detector_z[2 * i_central_mirror - i] = square_size_z;
+ l_aerogel_z[2 * i_central_mirror - i] = std::sqrt(1.0 + m_val * m_val) * R_min * square_size_z / (std::sqrt(1.0 + m_val * m_val) * R_max - m_val * square_size_z);
+ T_r_plus_g[2 * i_central_mirror - i] = std::sqrt(1.0 + m_val * m_val) * (R_max - R_min) - m_val / 2.0 * (square_size_z + l_aerogel_z[i]);
+ mProjectiveLengthInner = R_min * t; // <-- At the end of the loop this will be the maximum Z
+ }
+ // Coordinate radiali layer considerati
+ std::vector R0_detector;
+ std::vector R0_aerogel;
+ R0_detector.resize(number_of_sectors_in_z);
+ R0_aerogel.resize(number_of_sectors_in_z);
+ for (int i = 0; i < number_of_sectors_in_z; i++) {
+ R0_detector[i] = R0_tilt[i];
+ det_centers[i].SetXYZ(R0_detector[i], 0, R0_detector[i] * std::tan(theta_bi[i]));
+ R0_aerogel[i] = R0_tilt[i] - (T_r_plus_g[i]) * std::cos(theta_bi[i]);
+ rad_centers[i].SetXYZ(R0_aerogel[i], 0, R0_aerogel[i] * std::tan(theta_bi[i]));
+ angle_centers[i] = theta_bi[i];
+ }
+ }
+
+ void init(o2::framework::InitContext& initContext)
+ {
+ pRandomNumberGenerator.SetSeed(0); // fully randomize
+
+ // Load LUT for pt and eta smearing
+ if (flagIncludeTrackAngularRes && flagRICHLoadDelphesLUTs) {
+ std::map mapPdgLut;
+ const char* lutElChar = lutEl->c_str();
+ const char* lutMuChar = lutMu->c_str();
+ const char* lutPiChar = lutPi->c_str();
+ const char* lutKaChar = lutKa->c_str();
+ const char* lutPrChar = lutPr->c_str();
+
+ LOGF(info, "Will load electron lut file ..: %s for RICH PID", lutElChar);
+ LOGF(info, "Will load muon lut file ......: %s for RICH PID", lutMuChar);
+ LOGF(info, "Will load pion lut file ......: %s for RICH PID", lutPiChar);
+ LOGF(info, "Will load kaon lut file ......: %s for RICH PID", lutKaChar);
+ LOGF(info, "Will load proton lut file ....: %s for RICH PID", lutPrChar);
+
+ mapPdgLut.insert(std::make_pair(11, lutElChar));
+ mapPdgLut.insert(std::make_pair(13, lutMuChar));
+ mapPdgLut.insert(std::make_pair(211, lutPiChar));
+ mapPdgLut.insert(std::make_pair(321, lutKaChar));
+ mapPdgLut.insert(std::make_pair(2212, lutPrChar));
+
+ for (auto e : mapPdgLut) {
+ if (!mSmearer.loadTable(e.first, e.second)) {
+ LOG(fatal) << "Having issue with loading the LUT " << e.first << " " << e.second;
+ }
+ }
+ }
+
+ if (doQAplots) {
+ const AxisSpec axisMomentum{static_cast(nBinsP), 0.0f, +20.0f, "#it{p} (GeV/#it{c})"};
+ const AxisSpec axisAngle{static_cast(nBinsThetaRing), 0.0f, +0.30f, "Measured Cherenkov angle (rad)"};
+ histos.add("h2dAngleVsMomentumBarrelRICH", "h2dAngleVsMomentumBarrelRICH", kTH2F, {axisMomentum, axisAngle});
+
+ std::string particle_names1[5] = {"#it{e}", "#it{#mu}", "#it{#pi}", "#it{K}", "#it{p}"};
+ std::string particle_names2[5] = {"Elec", "Muon", "Pion", "Kaon", "Prot"};
+ for (int i_true = 0; i_true < 5; i_true++) {
+ std::string name_title_barrel_track_res = "h2dBarrelAngularResTrack" + particle_names2[i_true] + "VsP";
+ std::string name_title_barrel_total_res = "h2dBarrelAngularResTotal" + particle_names2[i_true] + "VsP";
+ const AxisSpec axisTrackAngularRes{static_cast(nBinsAngularRes), 0.0f, +5.0f, "Track angular resolution - " + particle_names1[i_true] + " (mrad)"};
+ const AxisSpec axisTotalAngularRes{static_cast(nBinsAngularRes), 0.0f, +5.0f, "Total angular resolution - " + particle_names1[i_true] + " (mrad)"};
+ histos.add(name_title_barrel_track_res.c_str(), name_title_barrel_track_res.c_str(), kTH2F, {axisMomentum, axisTrackAngularRes});
+ histos.add(name_title_barrel_total_res.c_str(), name_title_barrel_total_res.c_str(), kTH2F, {axisMomentum, axisTotalAngularRes});
+ }
+ for (int i_true = 0; i_true < 5; i_true++) {
+ for (int i_hyp = 0; i_hyp < 5; i_hyp++) {
+ std::string name_title = "h2dBarrelNsigmaTrue" + particle_names2[i_true] + "Vs" + particle_names2[i_hyp] + "Hypothesis";
+ if (i_true == i_hyp) {
+ const AxisSpec axisNsigmaCorrect{static_cast(nBinsNsigmaCorrectSpecies), -10.0f, +10.0f, "N#sigma - True " + particle_names1[i_true] + " vs " + particle_names1[i_hyp] + " hypothesis"};
+ histos.add(name_title.c_str(), name_title.c_str(), kTH2F, {axisMomentum, axisNsigmaCorrect});
+ } else {
+ const AxisSpec axisNsigmaWrong{static_cast(nBinsNsigmaWrongSpecies), -10.0f, +10.0f, "N#sigma - True " + particle_names1[i_true] + " vs " + particle_names1[i_hyp] + " hypothesis"};
+ histos.add(name_title.c_str(), name_title.c_str(), kTH2F, {axisMomentum, axisNsigmaWrong});
+ }
+ }
+ }
+ }
+
+ // Update projective parameters
+ updateProjectiveParameters();
+ }
+
+ /// Function to convert a McParticle into a perfect Track
+ /// \param particle the particle to convert (mcParticle)
+ /// \param o2track the address of the resulting TrackParCov
+ template
+ o2::track::TrackParCov convertMCParticleToO2Track(McParticleType& particle)
+ {
+ // FIXME: this is a fundamentally important piece of code.
+ // It could be placed in a utility file instead of here.
+ auto pdgInfo = pdg->GetParticle(particle.pdgCode());
+ int charge = 0;
+ if (pdgInfo != nullptr) {
+ charge = pdgInfo->Charge() / 3;
+ }
+ std::array params;
+ std::array covm = {0.};
+ float s, c, x;
+ o2::math_utils::sincos(particle.phi(), s, c);
+ o2::math_utils::rotateZInv(particle.vx(), particle.vy(), x, params[0], s, c);
+ params[1] = particle.vz();
+ params[2] = 0.; // since alpha = phi
+ auto theta = 2. * std::atan(std::exp(-particle.eta()));
+ params[3] = 1. / std::tan(theta);
+ params[4] = charge / particle.pt();
+
+ // Return TrackParCov
+ return o2::track::TrackParCov(x, particle.phi(), params, covm);
+ }
+
+ /// check if particle reaches radiator
+ /// \param track the input track
+ /// \param radius the radius of the layer you're calculating the length to
+ /// \param magneticField the magnetic field to use when propagating
+ bool checkMagfieldLimit(o2::track::TrackParCov track, float radius, float magneticField)
+ {
+ o2::math_utils::CircleXYf_t trcCircle;
+ float sna, csa;
+ track.getCircleParams(magneticField, trcCircle, sna, csa);
+
+ // distance between circle centers (one circle is at origin -> easy)
+ float centerDistance = std::hypot(trcCircle.xC, trcCircle.yC);
+
+ // condition of circles touching - if not satisfied returned value if false
+ if (centerDistance < trcCircle.rC + radius && centerDistance > std::fabs(trcCircle.rC - radius)) {
+ return true;
+ } else {
+ return false;
+ }
+ }
+
+ /// returns radiator radius in cm at considered eta (accounting for sector inclination)
+ /// \param eta the pseudorapidity of the tarck (assuming primary vertex at origin)
+ float radiusRipple(float eta)
+ {
+ float polar = 2.0 * std::atan(std::exp(-eta));
+ float i_sector = 0;
+ bool flag_sector = false;
+ for (int j_sec = 0; j_sec < mNumberSectors; j_sec++) {
+ if (polar > theta_max[j_sec] && polar < theta_min[j_sec]) {
+ flag_sector = true;
+ i_sector = j_sec;
+ }
+ }
+ if (flag_sector) {
+ float R_sec_rich = rad_centers[i_sector].X();
+ float z_sec_rich = rad_centers[i_sector].Z();
+ return (std::pow(R_sec_rich, 2) + std::pow(z_sec_rich, 2)) / (R_sec_rich + z_sec_rich / std::tan(polar));
+ } else {
+ return error_value;
+ }
+ }
+
+ /// returns Cherenkov angle in rad (above threshold) or bad flag (below threshold)
+ /// \param momentum the momentum of the tarck
+ /// \param mass the mass of the particle
+ float CherenkovAngle(float momentum, float mass)
+ {
+ // Check if particle is above the threshold
+ if (momentum > mass / std::sqrt(bRichRefractiveIndex * bRichRefractiveIndex - 1.0)) {
+ // Calculate angle
+ float angle = std::acos(std::sqrt(momentum * momentum + mass * mass) / (momentum * bRichRefractiveIndex));
+
+ // Mean number of detected photons (SiPM PDE is accountd in multiplicative factor!!!)
+ float meanNumberofDetectedPhotons = 230. * std::sin(angle) * std::sin(angle) * bRichRadiatorThickness;
+
+ // Require at least 3 photons on average for real angle reconstruction
+ if (meanNumberofDetectedPhotons > 3.) {
+ return angle;
+ } else {
+ return error_value;
+ }
+ } else {
+ return error_value;
+ }
+ }
+
+ /// returns linear interpolation
+ /// \param x the eta we want the resolution for
+ /// \param x0 the closest smaller available eta
+ /// \param x1 the closest larger available eta
+ /// \param y0 the resolution corresponding to x0
+ /// \param y1 the resolution corresponding to x1
+ float interpolate(double x, double x0, double x1, double y0, double y1)
+ {
+ float y = y0 + ((y1 - y0) / (x1 - x0)) * (x - x0);
+ return y;
+ }
+
+ /// returns angular resolution for considered track eta
+ /// \param eta the pseudorapidity of the tarck (assuming primary vertex at origin)
+ float AngularResolution(float eta)
+ {
+ // Vectors for sampling (USE ANALYTICAL EXTRAPOLATION FOR BETTER RESULTS)
+ float eta_sampling[] = {-2.000000, -1.909740, -1.731184, -1.552999, -1.375325, -1.198342, -1.022276, -0.847390, -0.673976, -0.502324, -0.332683, -0.165221, 0.000000, 0.165221, 0.332683, 0.502324, 0.673976, 0.847390, 1.022276, 1.198342, 1.375325, 1.552999, 1.731184, 1.909740, 2.000000};
+ float res_ring_sampling_with_abs_walls[] = {0.0009165, 0.000977, 0.001098, 0.001198, 0.001301, 0.001370, 0.001465, 0.001492, 0.001498, 0.001480, 0.001406, 0.001315, 0.001241, 0.001325, 0.001424, 0.001474, 0.001480, 0.001487, 0.001484, 0.001404, 0.001273, 0.001197, 0.001062, 0.000965, 0.0009165};
+ float res_ring_sampling_without_abs_walls[] = {0.0009165, 0.000977, 0.001095, 0.001198, 0.001300, 0.001369, 0.001468, 0.001523, 0.001501, 0.001426, 0.001299, 0.001167, 0.001092, 0.001179, 0.001308, 0.001407, 0.001491, 0.001508, 0.001488, 0.001404, 0.001273, 0.001196, 0.001061, 0.000965, 0.0009165};
+ int size_res_vector = sizeof(eta_sampling) / sizeof(eta_sampling[0]);
+ // Use binary search to find the lower and upper indices
+ int lowerIndex = std::lower_bound(eta_sampling, eta_sampling + size_res_vector, eta) - eta_sampling - 1;
+ int upperIndex = lowerIndex + 1;
+ if (lowerIndex >= 0 && upperIndex < size_res_vector) {
+ // Resolution interpolation
+ if (bRichFlagAbsorbingWalls) {
+ float interpolatedResRing = interpolate(eta, eta_sampling[lowerIndex], eta_sampling[upperIndex], res_ring_sampling_with_abs_walls[lowerIndex], res_ring_sampling_with_abs_walls[upperIndex]);
+ // std::cout << "Interpolated y value: " << interpolatedY << std::endl;
+ return interpolatedResRing;
+ } else {
+ float interpolatedResRing = interpolate(eta, eta_sampling[lowerIndex], eta_sampling[upperIndex], res_ring_sampling_without_abs_walls[lowerIndex], res_ring_sampling_without_abs_walls[upperIndex]);
+ // std::cout << "Interpolated y value: " << interpolatedY << std::endl;
+ return interpolatedResRing;
+ }
+ } else {
+ // std::cout << "Unable to interpolate. Target x value is outside the range of available data." << std::endl;
+ return error_value;
+ }
+ }
+
+ /// returns track angular resolution
+ /// \param pt the transverse momentum of the tarck
+ /// \param eta the pseudorapidity of the tarck
+ /// \param track_pt_resolution the absolute resolution on pt
+ /// \param track_pt_resolution the absolute resolution on eta
+ /// \param mass the mass of the particle
+ /// \param refractive_index the refractive index of the radiator
+ double calculate_track_time_resolution_advanced(float pt, float eta, float track_pt_resolution, float track_eta_resolution, float mass, float refractive_index)
+ {
+ // Compute tracking contribution to timing using the error propagation formula
+ // Uses light speed in m/ps, magnetic field in T (*0.1 for conversion kGauss -> T)
+ double a0 = mass * mass;
+ double a1 = refractive_index;
+ double dtheta_on_dpt = a0 / (pt * std::sqrt(a0 + std::pow(pt * std::cosh(eta), 2)) * std::sqrt(std::pow(pt * std::cosh(eta), 2) * (std::pow(a1, 2) - 1.0) - a0));
+ double dtheta_on_deta = (a0 * std::tanh(eta)) / (std::sqrt(a0 + std::pow(pt * std::cosh(eta), 2)) * std::sqrt(std::pow(pt * std::cosh(eta), 2) * (std::pow(a1, 2) - 1.0) - a0));
+ double track_angular_resolution = std::hypot(std::fabs(dtheta_on_dpt) * track_pt_resolution, std::fabs(dtheta_on_deta) * track_eta_resolution);
+ return track_angular_resolution;
+ }
+
+ void process(soa::Join::iterator const& collision, soa::Join const& tracks, aod::McParticles const&, aod::McCollisions const&)
+ {
+
+ o2::dataformats::VertexBase pvVtx({collision.posX(), collision.posY(), collision.posZ()},
+ {collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ()});
+
+ std::array mcPvCov = {0.};
+ o2::dataformats::VertexBase mcPvVtx({0.0f, 0.0f, 0.0f}, mcPvCov);
+ if (collision.has_mcCollision()) {
+ auto mcCollision = collision.mcCollision();
+ mcPvVtx.setX(mcCollision.posX());
+ mcPvVtx.setY(mcCollision.posY());
+ mcPvVtx.setZ(mcCollision.posZ());
+ } // else remains untreated for now
+
+ // First we compute the number of charged particles in the event
+ float dNdEta = 0.f;
+ if (flagRICHLoadDelphesLUTs) {
+ for (const auto& track : tracks) {
+ if (!track.has_mcParticle())
+ continue;
+ auto mcParticle = track.mcParticle();
+ if (std::abs(mcParticle.eta()) > multiplicityEtaRange) {
+ continue;
+ }
+ if (mcParticle.has_daughters()) {
+ continue;
+ }
+ const auto& pdgInfo = pdg->GetParticle(mcParticle.pdgCode());
+ if (!pdgInfo) {
+ // LOG(warning) << "PDG code " << mcParticle.pdgCode() << " not found in the database";
+ continue;
+ }
+ if (pdgInfo->Charge() == 0) {
+ continue;
+ }
+ dNdEta += 1.f;
+ }
+ }
+
+ for (const auto& track : tracks) {
+ // first step: find precise arrival time (if any)
+ // --- convert track into perfect track
+ if (!track.has_mcParticle()) // should always be OK but check please
+ continue;
+
+ auto mcParticle = track.mcParticle();
+ o2::track::TrackParCov o2track = convertMCParticleToO2Track(mcParticle);
+
+ // float xPv = error_value;
+ if (o2track.propagateToDCA(mcPvVtx, dBz)) {
+ // xPv = o2track.getX();
+ }
+
+ // get particle to calculate Cherenkov angle and resolution
+ auto pdgInfo = pdg->GetParticle(mcParticle.pdgCode());
+ if (pdgInfo == nullptr) {
+ continue;
+ }
+ float expectedAngleBarrelRich = CherenkovAngle(o2track.getP(), pdgInfo->Mass());
+ float barrelRICHAngularResolution = AngularResolution(o2track.getEta());
+ float projectiveRadiatorRadius = radiusRipple(o2track.getEta());
+ bool flagReachesRadiator = false;
+ if (projectiveRadiatorRadius > error_value + 1.) {
+ flagReachesRadiator = checkMagfieldLimit(o2track, projectiveRadiatorRadius, dBz);
+ }
+ /// DISCLAIMER: Exact extrapolation of angular resolution would require track propagation
+ /// to the RICH radiator (accounting sector inclination) in terms of (R,z).
+ /// The extrapolation with Eta is correct only if the primary vertex is at origin.
+ /// Discrepancies may be negligible, but would be more rigorous if propagation tool is available
+
+ // Smear with expected resolutions
+ float measuredAngleBarrelRich = pRandomNumberGenerator.Gaus(expectedAngleBarrelRich, barrelRICHAngularResolution);
+
+ // Now we calculate the expected arrival time following certain mass hypotheses
+ // and the (imperfect!) reconstructed track parametrizations
+ auto recoTrack = getTrackParCov(track);
+ if (recoTrack.propagateToDCA(pvVtx, dBz)) {
+ // xPv = recoTrack.getX();
+ }
+
+ // Straight to Nsigma
+ float deltaThetaBarrelRich[5], nSigmaBarrelRich[5];
+ int lpdg_array[5] = {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton};
+ float masses[5];
+
+ for (int ii = 0; ii < 5; ii++) {
+ nSigmaBarrelRich[ii] = error_value;
+
+ auto pdgInfoThis = pdg->GetParticle(lpdg_array[ii]);
+ masses[ii] = pdgInfoThis->Mass();
+ float hypothesisAngleBarrelRich = CherenkovAngle(recoTrack.getP(), masses[ii]);
+
+ // Evaluate total sigma (layer + tracking resolution)
+ float barrelTotalAngularReso = barrelRICHAngularResolution;
+ if (flagIncludeTrackAngularRes) {
+ double pt_resolution = std::pow(recoTrack.getP() / std::cosh(recoTrack.getEta()), 2) * std::sqrt(recoTrack.getSigma1Pt2());
+ double eta_resolution = std::fabs(std::sin(2.0 * std::atan(std::exp(-recoTrack.getEta())))) * std::sqrt(recoTrack.getSigmaTgl2());
+ if (flagRICHLoadDelphesLUTs) {
+ pt_resolution = mSmearer.getAbsPtRes(pdgInfoThis->PdgCode(), dNdEta, recoTrack.getEta(), recoTrack.getP() / std::cosh(recoTrack.getEta()));
+ eta_resolution = mSmearer.getAbsEtaRes(pdgInfoThis->PdgCode(), dNdEta, recoTrack.getEta(), recoTrack.getP() / std::cosh(recoTrack.getEta()));
+ }
+ // cout << endl << "Pt resolution: " << pt_resolution << ", Eta resolution: " << eta_resolution << endl << endl;
+ float barrelTrackAngularReso = calculate_track_time_resolution_advanced(recoTrack.getP() / std::cosh(recoTrack.getEta()), recoTrack.getEta(), pt_resolution, eta_resolution, masses[ii], bRichRefractiveIndex);
+ barrelTotalAngularReso = std::hypot(barrelRICHAngularResolution, barrelTrackAngularReso);
+ if (doQAplots && hypothesisAngleBarrelRich > error_value + 1. && measuredAngleBarrelRich > error_value + 1. && barrelRICHAngularResolution > error_value + 1. && flagReachesRadiator) {
+ float momentum = recoTrack.getP();
+ // float pseudorapidity = recoTrack.getEta();
+ // float transverse_momentum = momentum / std::cosh(pseudorapidity);
+ if (ii == 0 && std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[0])->PdgCode()) {
+ histos.fill(HIST("h2dBarrelAngularResTrackElecVsP"), momentum, 1000.0 * barrelTrackAngularReso);
+ histos.fill(HIST("h2dBarrelAngularResTotalElecVsP"), momentum, 1000.0 * barrelTotalAngularReso);
+ }
+ if (ii == 1 && std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[1])->PdgCode()) {
+ histos.fill(HIST("h2dBarrelAngularResTrackMuonVsP"), momentum, 1000.0 * barrelTrackAngularReso);
+ histos.fill(HIST("h2dBarrelAngularResTotalMuonVsP"), momentum, 1000.0 * barrelTotalAngularReso);
+ }
+ if (ii == 2 && std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[2])->PdgCode()) {
+ histos.fill(HIST("h2dBarrelAngularResTrackPionVsP"), momentum, 1000.0 * barrelTrackAngularReso);
+ histos.fill(HIST("h2dBarrelAngularResTotalPionVsP"), momentum, 1000.0 * barrelTotalAngularReso);
+ }
+ if (ii == 3 && std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[3])->PdgCode()) {
+ histos.fill(HIST("h2dBarrelAngularResTrackKaonVsP"), momentum, 1000.0 * barrelTrackAngularReso);
+ histos.fill(HIST("h2dBarrelAngularResTotalKaonVsP"), momentum, 1000.0 * barrelTotalAngularReso);
+ }
+ if (ii == 4 && std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[4])->PdgCode()) {
+ histos.fill(HIST("h2dBarrelAngularResTrackProtVsP"), momentum, 1000.0 * barrelTrackAngularReso);
+ histos.fill(HIST("h2dBarrelAngularResTotalProtVsP"), momentum, 1000.0 * barrelTotalAngularReso);
+ }
+ }
+ }
+
+ /// DISCLAIMER: here tracking is accounted only for momentum value, but not for track parameters at impact point on the
+ /// RICH radiator, since exact resolution would require photon generation and transport to photodetector.
+ /// Effects are expected to be negligible (a few tenths of a milliradian) but further studies are required !
+ if (hypothesisAngleBarrelRich > error_value + 1. && measuredAngleBarrelRich > error_value + 1. && barrelRICHAngularResolution > error_value + 1. && flagReachesRadiator) {
+ deltaThetaBarrelRich[ii] = hypothesisAngleBarrelRich - measuredAngleBarrelRich;
+ nSigmaBarrelRich[ii] = deltaThetaBarrelRich[ii] / barrelTotalAngularReso;
+ }
+ }
+
+ // Fill histograms
+ if (doQAplots) {
+ float momentum = recoTrack.getP();
+ float barrelRichTheta = measuredAngleBarrelRich;
+
+ if (barrelRichTheta > error_value + 1. && barrelRICHAngularResolution > error_value + 1. && flagReachesRadiator) {
+ histos.fill(HIST("h2dAngleVsMomentumBarrelRICH"), momentum, barrelRichTheta);
+
+ if (std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[0])->PdgCode()) {
+ histos.fill(HIST("h2dBarrelNsigmaTrueElecVsElecHypothesis"), momentum, nSigmaBarrelRich[0]);
+ histos.fill(HIST("h2dBarrelNsigmaTrueElecVsMuonHypothesis"), momentum, nSigmaBarrelRich[1]);
+ histos.fill(HIST("h2dBarrelNsigmaTrueElecVsPionHypothesis"), momentum, nSigmaBarrelRich[2]);
+ histos.fill(HIST("h2dBarrelNsigmaTrueElecVsKaonHypothesis"), momentum, nSigmaBarrelRich[3]);
+ histos.fill(HIST("h2dBarrelNsigmaTrueElecVsProtHypothesis"), momentum, nSigmaBarrelRich[4]);
+ }
+ if (std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[1])->PdgCode()) {
+ histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsElecHypothesis"), momentum, nSigmaBarrelRich[0]);
+ histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsMuonHypothesis"), momentum, nSigmaBarrelRich[1]);
+ histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsPionHypothesis"), momentum, nSigmaBarrelRich[2]);
+ histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsKaonHypothesis"), momentum, nSigmaBarrelRich[3]);
+ histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsProtHypothesis"), momentum, nSigmaBarrelRich[4]);
+ }
+ if (std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[2])->PdgCode()) {
+ histos.fill(HIST("h2dBarrelNsigmaTruePionVsElecHypothesis"), momentum, nSigmaBarrelRich[0]);
+ histos.fill(HIST("h2dBarrelNsigmaTruePionVsMuonHypothesis"), momentum, nSigmaBarrelRich[1]);
+ histos.fill(HIST("h2dBarrelNsigmaTruePionVsPionHypothesis"), momentum, nSigmaBarrelRich[2]);
+ histos.fill(HIST("h2dBarrelNsigmaTruePionVsKaonHypothesis"), momentum, nSigmaBarrelRich[3]);
+ histos.fill(HIST("h2dBarrelNsigmaTruePionVsProtHypothesis"), momentum, nSigmaBarrelRich[4]);
+ }
+ if (std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[3])->PdgCode()) {
+ histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsElecHypothesis"), momentum, nSigmaBarrelRich[0]);
+ histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsMuonHypothesis"), momentum, nSigmaBarrelRich[1]);
+ histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsPionHypothesis"), momentum, nSigmaBarrelRich[2]);
+ histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsKaonHypothesis"), momentum, nSigmaBarrelRich[3]);
+ histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsProtHypothesis"), momentum, nSigmaBarrelRich[4]);
+ }
+ if (std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[4])->PdgCode()) {
+ histos.fill(HIST("h2dBarrelNsigmaTrueProtVsElecHypothesis"), momentum, nSigmaBarrelRich[0]);
+ histos.fill(HIST("h2dBarrelNsigmaTrueProtVsMuonHypothesis"), momentum, nSigmaBarrelRich[1]);
+ histos.fill(HIST("h2dBarrelNsigmaTrueProtVsPionHypothesis"), momentum, nSigmaBarrelRich[2]);
+ histos.fill(HIST("h2dBarrelNsigmaTrueProtVsKaonHypothesis"), momentum, nSigmaBarrelRich[3]);
+ histos.fill(HIST("h2dBarrelNsigmaTrueProtVsProtHypothesis"), momentum, nSigmaBarrelRich[4]);
+ }
+ }
+ }
+
+ // Sigmas have been fully calculated. Please populate the NSigma helper table (once per track)
+ upgradeRich(nSigmaBarrelRich[0], nSigmaBarrelRich[1], nSigmaBarrelRich[2], nSigmaBarrelRich[3], nSigmaBarrelRich[4]);
+ }
+ }
+};
+
+WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
+{
+ return WorkflowSpec{adaptAnalysisTask(cfgc)};
+}
diff --git a/ALICE3/TableProducer/OTF/onTheFlyTOFPID.cxx b/ALICE3/TableProducer/OTF/onTheFlyTOFPID.cxx
new file mode 100644
index 00000000000..f9d0aa4cb7e
--- /dev/null
+++ b/ALICE3/TableProducer/OTF/onTheFlyTOFPID.cxx
@@ -0,0 +1,604 @@
+// 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.
+
+//
+// Task to add a table of track parameters propagated to the primary vertex
+//
+
+#include
+#include "Framework/AnalysisDataModel.h"
+#include "Framework/AnalysisTask.h"
+#include "Framework/runDataProcessing.h"
+#include "Framework/RunningWorkflowInfo.h"
+#include "Framework/HistogramRegistry.h"
+#include "Framework/O2DatabasePDGPlugin.h"
+#include "Framework/ASoAHelpers.h"
+#include "Common/DataModel/TrackSelectionTables.h"
+#include "Common/Core/trackUtilities.h"
+#include "ReconstructionDataFormats/DCA.h"
+#include "DetectorsBase/Propagator.h"
+#include "DetectorsBase/GeometryManager.h"
+#include "CommonUtils/NameConf.h"
+#include "CCDB/CcdbApi.h"
+#include "CCDB/BasicCCDBManager.h"
+#include "DataFormatsParameters/GRPMagField.h"
+#include "DataFormatsCalibration/MeanVertexObject.h"
+#include "CommonConstants/GeomConstants.h"
+#include "CommonConstants/PhysicsConstants.h"
+#include "TRandom3.h"
+#include "ALICE3/DataModel/OTFTOF.h"
+#include "DetectorsVertexing/HelixHelper.h"
+#include "TableHelper.h"
+#include "ALICE3/Core/DelphesO2TrackSmearer.h"
+
+/// \file onTheFlyTOFPID.cxx
+///
+/// \brief This task goes straight from a combination of track table and mcParticles
+/// and a custom TOF configuration to a table of TOF NSigmas for the particles
+/// being analysed. It currently contemplates 5 particle types:
+/// electrons, pions, kaons, protons and muons
+///
+/// More particles could be added but would have to be added to the LUT
+/// being used in the onTheFly tracker task.
+///
+/// \author David Dobrigkeit Chinellato, UNICAMP, Nicola Nicassio, University and INFN Bari
+
+using namespace o2;
+using namespace o2::framework;
+
+struct OnTheFlyTOFPID {
+ Produces upgradeTof;
+
+ // necessary for particle charges
+ Service pdg;
+
+ // these are the settings governing the TOF layers to be used
+ // note that there are two layers foreseen for now: inner and outer TOF
+ // more could be added (especially a disk TOF at a certain z?)
+ // in the evolution of this effort
+ Configurable dBz{"dBz", 20, "magnetic field (kilogauss)"};
+ Configurable innerTOFRadius{"innerTOFRadius", 20, "barrel inner TOF radius (cm)"};
+ Configurable outerTOFRadius{"outerTOFRadius", 80, "barrel outer TOF radius (cm)"};
+ Configurable innerTOFTimeReso{"innerTOFTimeReso", 20, "barrel inner TOF time error (ps)"};
+ Configurable outerTOFTimeReso{"outerTOFTimeReso", 20, "barrel outer TOF time error (ps)"};
+ Configurable nStepsLIntegrator{"nStepsLIntegrator", 200, "number of steps in length integrator"};
+ Configurable doQAplots{"doQAplots", true, "do basic velocity plot qa"};
+ Configurable nBinsBeta{"nBinsBeta", 2200, "number of bins in beta"};
+ Configurable nBinsP{"nBinsP", 80, "number of bins in momentum"};
+ Configurable nBinsTrackLengthInner{"nBinsTrackLengthInner", 300, "number of bins in track length"};
+ Configurable nBinsTrackLengthOuter{"nBinsTrackLengthOuter", 300, "number of bins in track length"};
+ Configurable nBinsTrackDeltaLength{"nBinsTrackDeltaLength", 100, "number of bins in delta track length"};
+ Configurable nBinsNsigmaCorrectSpecies{"nBinsNsigmaCorrectSpecies", 200, "number of bins in Nsigma plot (correct speies)"};
+ Configurable nBinsNsigmaWrongSpecies{"nBinsNsigmaWrongSpecies", 200, "number of bins in Nsigma plot (wrong species)"};
+ Configurable nBinsTimeRes{"nBinsTimeRes", 400, "number of bins plots time resolution"};
+ Configurable nBinsRelativeEtaPt{"nBinsRelativeEtaPt", 400, "number of bins plots pt and eta relative errors"};
+ Configurable nBinsEta{"nBinsEta", 400, "number of bins plot relative eta error"};
+ Configurable flagIncludeTrackTimeRes{"flagIncludeTrackTimeRes", true, "flag to include or exclude track time resolution"};
+ Configurable multiplicityEtaRange{"multiplicityEtaRange", 0.800000012, "eta range to compute the multiplicity"};
+ Configurable flagTOFLoadDelphesLUTs{"flagTOFLoadDelphesLUTs", false, "flag to load Delphes LUTs for tracking correction (use recoTrack parameters if false)"};
+
+ Configurable lutEl{"lutEl", "lutCovm.el.dat", "LUT for electrons"};
+ Configurable lutMu{"lutMu", "lutCovm.mu.dat", "LUT for muons"};
+ Configurable lutPi{"lutPi", "lutCovm.pi.dat", "LUT for pions"};
+ Configurable lutKa{"lutKa", "lutCovm.ka.dat", "LUT for kaons"};
+ Configurable lutPr{"lutPr", "lutCovm.pr.dat", "LUT for protons"};
+
+ o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrNONE;
+
+ // Track smearer (here used to get absolute pt and eta uncertainties if flagTOFLoadDelphesLUTs is true)
+ o2::delphes::DelphesO2TrackSmearer mSmearer;
+
+ // needed: random number generator for smearing
+ TRandom3 pRandomNumberGenerator;
+
+ // for handling basic QA histograms if requested
+ HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject};
+
+ void init(o2::framework::InitContext& initContext)
+ {
+ pRandomNumberGenerator.SetSeed(0); // fully randomize
+
+ // Load LUT for pt and eta smearing
+ if (flagIncludeTrackTimeRes && flagTOFLoadDelphesLUTs) {
+ std::map mapPdgLut;
+ const char* lutElChar = lutEl->c_str();
+ const char* lutMuChar = lutMu->c_str();
+ const char* lutPiChar = lutPi->c_str();
+ const char* lutKaChar = lutKa->c_str();
+ const char* lutPrChar = lutPr->c_str();
+
+ LOGF(info, "Will load electron lut file ..: %s for TOF PID", lutElChar);
+ LOGF(info, "Will load muon lut file ......: %s for TOF PID", lutMuChar);
+ LOGF(info, "Will load pion lut file ......: %s for TOF PID", lutPiChar);
+ LOGF(info, "Will load kaon lut file ......: %s for TOF PID", lutKaChar);
+ LOGF(info, "Will load proton lut file ....: %s for TOF PID", lutPrChar);
+
+ mapPdgLut.insert(std::make_pair(11, lutElChar));
+ mapPdgLut.insert(std::make_pair(13, lutMuChar));
+ mapPdgLut.insert(std::make_pair(211, lutPiChar));
+ mapPdgLut.insert(std::make_pair(321, lutKaChar));
+ mapPdgLut.insert(std::make_pair(2212, lutPrChar));
+
+ for (auto e : mapPdgLut) {
+ if (!mSmearer.loadTable(e.first, e.second)) {
+ LOG(fatal) << "Having issue with loading the LUT " << e.first << " " << e.second;
+ }
+ }
+ }
+
+ if (doQAplots) {
+ const AxisSpec axisMomentum{static_cast(nBinsP), 0.0f, +4.0f, "#it{p} (GeV/#it{c})"};
+ const AxisSpec axisMomentumSmall{static_cast(nBinsP), 0.0f, +1.0f, "#it{p} (GeV/#it{c})"};
+ const AxisSpec axisVelocity{static_cast(nBinsBeta), 0.0f, +1.1f, "Measured #beta"};
+ const AxisSpec axisTrackLengthInner{static_cast(nBinsTrackLengthInner), 0.0f, 60.0f, "Track length (cm)"};
+ const AxisSpec axisTrackLengthOuter{static_cast(nBinsTrackLengthOuter), 0.0f, 300.0f, "Track length (cm)"};
+ const AxisSpec axisTrackDeltaLength{static_cast(nBinsTrackDeltaLength), 0.0f, 30.0f, "Delta Track length (cm)"};
+ histos.add("h2dVelocityVsMomentumInner", "h2dVelocityVsMomentumInner", kTH2F, {axisMomentum, axisVelocity});
+ histos.add("h2dVelocityVsMomentumOuter", "h2dVelocityVsMomentumOuter", kTH2F, {axisMomentum, axisVelocity});
+ histos.add("h2dTrackLengthInnerVsPt", "h2dTrackLengthInnerVsPt", kTH2F, {axisMomentumSmall, axisTrackLengthInner});
+ histos.add("h2dTrackLengthOuterVsPt", "h2dTrackLengthOuterVsPt", kTH2F, {axisMomentumSmall, axisTrackLengthOuter});
+
+ histos.add("h2dTrackLengthInnerRecoVsPt", "h2dTrackLengthInnerRecoVsPt", kTH2F, {axisMomentumSmall, axisTrackLengthInner});
+ histos.add("h2dTrackLengthOuterRecoVsPt", "h2dTrackLengthOuterRecoVsPt", kTH2F, {axisMomentumSmall, axisTrackLengthOuter});
+
+ histos.add("h2dDeltaTrackLengthInnerVsPt", "h2dDeltaTrackLengthInnerVsPt", kTH2F, {axisMomentumSmall, axisTrackDeltaLength});
+ histos.add("h2dDeltaTrackLengthOuterVsPt", "h2dDeltaTrackLengthOuterVsPt", kTH2F, {axisMomentumSmall, axisTrackDeltaLength});
+
+ const AxisSpec axisPt{static_cast(nBinsP), 0.0f, +4.0f, "#it{p_{T}} (GeV/#it{c})"};
+ const AxisSpec axisEta{static_cast(nBinsEta), -2.0f, +2.0f, "#eta"};
+ const AxisSpec axisRelativePt{static_cast(nBinsRelativeEtaPt), 0.0f, +10.0f, "#it{#sigma_{p_{T}}} / #it{p_{T}} (%)"};
+ const AxisSpec axisRelativeEta{static_cast(nBinsRelativeEtaPt), 0.0f, +10.0f, "#it{#sigma_{#eta}} / #it{#eta} (%)"};
+ histos.add("h2dRelativePtResolution", "h2dRelativePtResolution", kTH2F, {axisPt, axisRelativePt});
+ histos.add("h2dRelativeEtaResolution", "h2dRelativeEtaResolution", kTH2F, {axisEta, axisRelativeEta});
+
+ std::string particle_names1[5] = {"#it{e}", "#it{#mu}", "#it{#pi}", "#it{K}", "#it{p}"};
+ std::string particle_names2[5] = {"Elec", "Muon", "Pion", "Kaon", "Prot"};
+ for (int i_true = 0; i_true < 5; i_true++) {
+ std::string name_title_inner_track_res = "h2dInnerTimeResTrack" + particle_names2[i_true] + "VsP";
+ std::string name_title_inner_total_res = "h2dInnerTimeResTotal" + particle_names2[i_true] + "VsP";
+ std::string name_title_outer_track_res = "h2dOuterTimeResTrack" + particle_names2[i_true] + "VsP";
+ std::string name_title_outer_total_res = "h2dOuterTimeResTotal" + particle_names2[i_true] + "VsP";
+ const AxisSpec axisTrackTimeRes{static_cast(nBinsTimeRes), 0.0f, +200.0f, "Track time resolution - " + particle_names1[i_true] + " (ps)"};
+ const AxisSpec axisTotalTimeRes{static_cast(nBinsTimeRes), 0.0f, +200.0f, "Total time resolution - " + particle_names1[i_true] + " (ps)"};
+ histos.add(name_title_inner_track_res.c_str(), name_title_inner_track_res.c_str(), kTH2F, {axisMomentum, axisTrackTimeRes});
+ histos.add(name_title_inner_total_res.c_str(), name_title_inner_total_res.c_str(), kTH2F, {axisMomentum, axisTotalTimeRes});
+ histos.add(name_title_outer_track_res.c_str(), name_title_outer_track_res.c_str(), kTH2F, {axisMomentum, axisTrackTimeRes});
+ histos.add(name_title_outer_total_res.c_str(), name_title_outer_total_res.c_str(), kTH2F, {axisMomentum, axisTotalTimeRes});
+ }
+
+ for (int i_true = 0; i_true < 5; i_true++) {
+ for (int i_hyp = 0; i_hyp < 5; i_hyp++) {
+ std::string name_title_inner = "h2dInnerNsigmaTrue" + particle_names2[i_true] + "Vs" + particle_names2[i_hyp] + "Hypothesis";
+ std::string name_title_outer = "h2dOuterNsigmaTrue" + particle_names2[i_true] + "Vs" + particle_names2[i_hyp] + "Hypothesis";
+ if (i_true == i_hyp) {
+ const AxisSpec axisNsigmaCorrect{static_cast(nBinsNsigmaCorrectSpecies), -10.0f, +10.0f, "N#sigma - True " + particle_names1[i_true] + " vs " + particle_names1[i_hyp] + " hypothesis"};
+ histos.add(name_title_inner.c_str(), name_title_inner.c_str(), kTH2F, {axisMomentum, axisNsigmaCorrect});
+ histos.add(name_title_outer.c_str(), name_title_outer.c_str(), kTH2F, {axisMomentum, axisNsigmaCorrect});
+ } else {
+ const AxisSpec axisNsigmaWrong{static_cast(nBinsNsigmaWrongSpecies), -10.0f, +10.0f, "N#sigma - True " + particle_names1[i_true] + " vs " + particle_names1[i_hyp] + " hypothesis"};
+ histos.add(name_title_inner.c_str(), name_title_inner.c_str(), kTH2F, {axisMomentum, axisNsigmaWrong});
+ histos.add(name_title_outer.c_str(), name_title_outer.c_str(), kTH2F, {axisMomentum, axisNsigmaWrong});
+ }
+ }
+ }
+ }
+ }
+
+ /// Function to convert a McParticle into a perfect Track
+ /// \param particle the particle to convert (mcParticle)
+ /// \param o2track the address of the resulting TrackParCov
+ template
+ o2::track::TrackParCov convertMCParticleToO2Track(McParticleType& particle)
+ {
+ // FIXME: this is a fundamentally important piece of code.
+ // It could be placed in a utility file instead of here.
+ auto pdgInfo = pdg->GetParticle(particle.pdgCode());
+ int charge = 0;
+ if (pdgInfo != nullptr) {
+ charge = pdgInfo->Charge() / 3;
+ }
+ std::array params;
+ std::array covm = {0.};
+ float s, c, x;
+ o2::math_utils::sincos(particle.phi(), s, c);
+ o2::math_utils::rotateZInv(particle.vx(), particle.vy(), x, params[0], s, c);
+ params[1] = particle.vz();
+ params[2] = 0.; // since alpha = phi
+ auto theta = 2. * std::atan(std::exp(-particle.eta()));
+ params[3] = 1. / std::tan(theta);
+ params[4] = charge / particle.pt();
+
+ // Return TrackParCov
+ return o2::track::TrackParCov(x, particle.phi(), params, covm);
+ }
+
+ /// function to calculate track length of this track up to a certain radius
+ /// \param track the input track
+ /// \param radius the radius of the layer you're calculating the length to
+ /// \param magneticField the magnetic field to use when propagating
+ float trackLength(o2::track::TrackParCov track, float radius, float magneticField)
+ {
+ // don't make use of the track parametrization
+ float length = -100;
+
+ o2::math_utils::CircleXYf_t trcCircle;
+ float sna, csa;
+ track.getCircleParams(magneticField, trcCircle, sna, csa);
+
+ // distance between circle centers (one circle is at origin -> easy)
+ float centerDistance = std::hypot(trcCircle.xC, trcCircle.yC);
+
+ // condition of circles touching - if not satisfied returned length will be -100
+ if (centerDistance < trcCircle.rC + radius && centerDistance > fabs(trcCircle.rC - radius)) {
+ length = 0.0f;
+
+ // base radical direction
+ float ux = trcCircle.xC / centerDistance;
+ float uy = trcCircle.yC / centerDistance;
+ // calculate perpendicular vector (normalized) for +/- displacement
+ float vx = -uy;
+ float vy = +ux;
+ // calculate coordinate for radical line
+ float radical = (centerDistance * centerDistance - trcCircle.rC * trcCircle.rC + radius * radius) / (2.0f * centerDistance);
+ // calculate absolute displacement from center-to-center axis
+ float displace = (0.5f / centerDistance) * TMath::Sqrt(
+ (-centerDistance + trcCircle.rC - radius) *
+ (-centerDistance - trcCircle.rC + radius) *
+ (-centerDistance + trcCircle.rC + radius) *
+ (centerDistance + trcCircle.rC + radius));
+
+ // possible intercept points of track and TOF layer in 2D plane
+ float point1[2] = {radical * ux + displace * vx, radical * uy + displace * vy};
+ float point2[2] = {radical * ux - displace * vx, radical * uy - displace * vy};
+
+ // decide on correct intercept point
+ std::array mom;
+ track.getPxPyPzGlo(mom);
+ float scalarProduct1 = point1[0] * mom[0] + point1[1] * mom[1];
+ float scalarProduct2 = point2[0] * mom[0] + point2[1] * mom[1];
+
+ // get start point
+ std::array startPoint;
+ track.getXYZGlo(startPoint);
+
+ float cosAngle = -1000, modulus = -1000;
+
+ if (scalarProduct1 > scalarProduct2) {
+ modulus = std::hypot(point1[0] - trcCircle.xC, point1[1] - trcCircle.yC) * std::hypot(startPoint[0] - trcCircle.xC, startPoint[1] - trcCircle.yC);
+ cosAngle = (point1[0] - trcCircle.xC) * (startPoint[0] - trcCircle.xC) + (point1[1] - trcCircle.yC) * (startPoint[0] - trcCircle.yC);
+ } else {
+ modulus = std::hypot(point2[0] - trcCircle.xC, point2[1] - trcCircle.yC) * std::hypot(startPoint[0] - trcCircle.xC, startPoint[1] - trcCircle.yC);
+ cosAngle = (point2[0] - trcCircle.xC) * (startPoint[0] - trcCircle.xC) + (point2[1] - trcCircle.yC) * (startPoint[0] - trcCircle.yC);
+ }
+ cosAngle /= modulus;
+ length = trcCircle.rC * TMath::ACos(cosAngle);
+ length *= sqrt(1.0f + track.getTgl() * track.getTgl());
+ }
+ return length;
+ }
+
+ /// returns velocity in centimeters per picoseconds
+ /// \param momentum the momentum of the tarck
+ /// \param mass the mass of the particle
+ float velocity(float momentum, float mass)
+ {
+ float a = std::pow(momentum / mass, 2);
+ // uses light speed in cm/ps so output is in those units
+ return (o2::constants::physics::LightSpeedCm2NS / 1e+3) * std::sqrt(a / (1 + a));
+ }
+
+ /// returns track time resolution
+ /// \param pt the transverse momentum of the tarck
+ /// \param eta the pseudorapidity of the tarck
+ /// \param track_pt_resolution the absolute resolution on pt
+ /// \param track_pt_resolution the absolute resolution on eta
+ /// \param mass the mass of the particle
+ /// \param det_radius the radius of the cylindrical layer
+ /// \param magneticField the magnetic field (along Z)
+ double calculate_track_time_resolution_advanced(float pt, float eta, float track_pt_resolution, float track_eta_resolution, float mass, float det_radius, float magneticField)
+ {
+ // Compute tracking contribution to timing using the error propagation formula
+ // Uses light speed in m/ps, magnetic field in T (*0.1 for conversion kGauss -> T)
+ double a0 = mass * mass;
+ double a1 = 0.299792458 * (0.1 * magneticField) * (0.01 * o2::constants::physics::LightSpeedCm2NS / 1e+3);
+ double a2 = (det_radius * 0.01) * (det_radius * 0.01) * (0.299792458) * (0.299792458) * (0.1 * magneticField) * (0.1 * magneticField) / 2.0;
+ double dtof_on_dpt = (std::pow(pt, 4) * std::pow(std::cosh(eta), 2) * std::acos(1.0 - a2 / std::pow(pt, 2)) - 2.0 * a2 * std::pow(pt, 2) * (a0 + std::pow(pt * std::cosh(eta), 2)) / std::sqrt(a2 * (2.0 * std::pow(pt, 2) - a2))) / (a1 * std::pow(pt, 3) * std::sqrt(a0 + std::pow(pt * std::cosh(eta), 2)));
+ double dtof_on_deta = std::pow(pt, 2) * std::sinh(eta) * std::cosh(eta) * std::acos(1.0 - a2 / std::pow(pt, 2)) / (a1 * std::sqrt(a0 + std::pow(pt * std::cosh(eta), 2)));
+ double track_time_resolution = std::hypot(std::fabs(dtof_on_dpt) * track_pt_resolution, std::fabs(dtof_on_deta) * track_eta_resolution);
+ return track_time_resolution;
+ }
+
+ void process(soa::Join::iterator const& collision, soa::Join const& tracks, aod::McParticles const&, aod::McCollisions const&)
+ {
+ o2::dataformats::VertexBase pvVtx({collision.posX(), collision.posY(), collision.posZ()},
+ {collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ()});
+
+ std::array mcPvCov = {0.};
+ o2::dataformats::VertexBase mcPvVtx({0.0f, 0.0f, 0.0f}, mcPvCov);
+ if (collision.has_mcCollision()) {
+ auto mcCollision = collision.mcCollision();
+ mcPvVtx.setX(mcCollision.posX());
+ mcPvVtx.setY(mcCollision.posY());
+ mcPvVtx.setZ(mcCollision.posZ());
+ } // else remains untreated for now
+
+ // First we compute the number of charged particles in the event if LUTs are loaded
+ float dNdEta = 0.f;
+ if (flagTOFLoadDelphesLUTs) {
+ for (const auto& track : tracks) {
+ if (!track.has_mcParticle())
+ continue;
+ auto mcParticle = track.mcParticle();
+ if (std::abs(mcParticle.eta()) > multiplicityEtaRange) {
+ continue;
+ }
+ if (mcParticle.has_daughters()) {
+ continue;
+ }
+ const auto& pdgInfo = pdg->GetParticle(mcParticle.pdgCode());
+ if (!pdgInfo) {
+ // LOG(warning) << "PDG code " << mcParticle.pdgCode() << " not found in the database";
+ continue;
+ }
+ if (pdgInfo->Charge() == 0) {
+ continue;
+ }
+ dNdEta += 1.f;
+ }
+ }
+
+ for (const auto& track : tracks) {
+ // first step: find precise arrival time (if any)
+ // --- convert track into perfect track
+ if (!track.has_mcParticle()) // should always be OK but check please
+ continue;
+
+ auto mcParticle = track.mcParticle();
+ o2::track::TrackParCov o2track = convertMCParticleToO2Track(mcParticle);
+
+ float xPv = -100, trackLengthInnerTOF = -1, trackLengthOuterTOF = -1;
+ if (o2track.propagateToDCA(mcPvVtx, dBz))
+ xPv = o2track.getX();
+ if (xPv > -99.) {
+ trackLengthInnerTOF = trackLength(o2track, innerTOFRadius, dBz);
+ trackLengthOuterTOF = trackLength(o2track, outerTOFRadius, dBz);
+ }
+
+ // get mass to calculate velocity
+ auto pdgInfo = pdg->GetParticle(mcParticle.pdgCode());
+ if (pdgInfo == nullptr) {
+ continue;
+ }
+ float expectedTimeInnerTOF = trackLengthInnerTOF / velocity(o2track.getP(), pdgInfo->Mass());
+ float expectedTimeOuterTOF = trackLengthOuterTOF / velocity(o2track.getP(), pdgInfo->Mass());
+
+ // Smear with expected resolutions
+ float measuredTimeInnerTOF = pRandomNumberGenerator.Gaus(expectedTimeInnerTOF, innerTOFTimeReso);
+ float measuredTimeOuterTOF = pRandomNumberGenerator.Gaus(expectedTimeOuterTOF, outerTOFTimeReso);
+
+ // Now we calculate the expected arrival time following certain mass hypotheses
+ // and the (imperfect!) reconstructed track parametrizations
+ float trackLengthRecoInnerTOF = -1, trackLengthRecoOuterTOF = -1;
+ auto recoTrack = getTrackParCov(track);
+ if (recoTrack.propagateToDCA(pvVtx, dBz))
+ xPv = recoTrack.getX();
+ if (xPv > -99.) {
+ trackLengthRecoInnerTOF = trackLength(recoTrack, innerTOFRadius, dBz);
+ trackLengthRecoOuterTOF = trackLength(recoTrack, outerTOFRadius, dBz);
+ }
+
+ // Straight to Nsigma
+ float deltaTimeInnerTOF[5], nSigmaInnerTOF[5];
+ float deltaTimeOuterTOF[5], nSigmaOuterTOF[5];
+ int lpdg_array[5] = {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton};
+ float masses[5];
+
+ if (doQAplots) {
+ float momentum = recoTrack.getP();
+ // unit conversion: length in cm, time in ps
+ float innerBeta = 1e+3 * (trackLengthInnerTOF / measuredTimeInnerTOF) / o2::constants::physics::LightSpeedCm2NS;
+ float outerBeta = 1e+3 * (trackLengthOuterTOF / measuredTimeOuterTOF) / o2::constants::physics::LightSpeedCm2NS;
+ if (trackLengthRecoInnerTOF > 0) {
+ histos.fill(HIST("h2dVelocityVsMomentumInner"), momentum, innerBeta);
+ histos.fill(HIST("h2dTrackLengthInnerVsPt"), o2track.getPt(), trackLengthInnerTOF);
+ histos.fill(HIST("h2dTrackLengthInnerRecoVsPt"), o2track.getPt(), trackLengthRecoInnerTOF);
+ }
+ if (trackLengthRecoOuterTOF > 0) {
+ histos.fill(HIST("h2dVelocityVsMomentumOuter"), momentum, outerBeta);
+ histos.fill(HIST("h2dTrackLengthOuterVsPt"), o2track.getPt(), trackLengthOuterTOF);
+ histos.fill(HIST("h2dTrackLengthOuterRecoVsPt"), o2track.getPt(), trackLengthRecoOuterTOF);
+ }
+ }
+
+ for (int ii = 0; ii < 5; ii++) {
+ nSigmaInnerTOF[ii] = -100;
+ nSigmaOuterTOF[ii] = -100;
+
+ auto pdgInfoThis = pdg->GetParticle(lpdg_array[ii]);
+ masses[ii] = pdgInfoThis->Mass();
+ deltaTimeInnerTOF[ii] = trackLengthRecoInnerTOF / velocity(recoTrack.getP(), masses[ii]) - measuredTimeInnerTOF;
+ deltaTimeOuterTOF[ii] = trackLengthRecoOuterTOF / velocity(recoTrack.getP(), masses[ii]) - measuredTimeOuterTOF;
+
+ // Evaluate total sigma (layer + tracking resolution)
+ float innerTotalTimeReso = innerTOFTimeReso;
+ float outerTotalTimeReso = outerTOFTimeReso;
+ if (flagIncludeTrackTimeRes) {
+ double pt_resolution = std::pow(recoTrack.getP() / std::cosh(recoTrack.getEta()), 2) * std::sqrt(recoTrack.getSigma1Pt2());
+ double eta_resolution = std::fabs(std::sin(2.0 * std::atan(std::exp(-recoTrack.getEta())))) * std::sqrt(recoTrack.getSigmaTgl2());
+ if (flagTOFLoadDelphesLUTs) {
+ pt_resolution = mSmearer.getAbsPtRes(pdgInfoThis->PdgCode(), dNdEta, recoTrack.getEta(), recoTrack.getP() / std::cosh(recoTrack.getEta()));
+ eta_resolution = mSmearer.getAbsEtaRes(pdgInfoThis->PdgCode(), dNdEta, recoTrack.getEta(), recoTrack.getP() / std::cosh(recoTrack.getEta()));
+ }
+ float innerTrackTimeReso = calculate_track_time_resolution_advanced(recoTrack.getP() / std::cosh(recoTrack.getEta()), recoTrack.getEta(), pt_resolution, eta_resolution, masses[ii], innerTOFRadius, dBz);
+ float outerTrackTimeReso = calculate_track_time_resolution_advanced(recoTrack.getP() / std::cosh(recoTrack.getEta()), recoTrack.getEta(), pt_resolution, eta_resolution, masses[ii], outerTOFRadius, dBz);
+ innerTotalTimeReso = std::hypot(innerTOFTimeReso, innerTrackTimeReso);
+ outerTotalTimeReso = std::hypot(outerTOFTimeReso, outerTrackTimeReso);
+
+ if (doQAplots && trackLengthRecoInnerTOF > 0) {
+ float momentum = recoTrack.getP();
+ if (ii == 0 && std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[0])->PdgCode()) {
+ histos.fill(HIST("h2dInnerTimeResTrackElecVsP"), momentum, innerTrackTimeReso);
+ histos.fill(HIST("h2dInnerTimeResTotalElecVsP"), momentum, innerTotalTimeReso);
+ }
+ if (ii == 1 && std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[1])->PdgCode()) {
+ histos.fill(HIST("h2dInnerTimeResTrackMuonVsP"), momentum, innerTrackTimeReso);
+ histos.fill(HIST("h2dInnerTimeResTotalMuonVsP"), momentum, innerTotalTimeReso);
+ }
+ if (ii == 2 && std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[2])->PdgCode()) {
+ histos.fill(HIST("h2dInnerTimeResTrackPionVsP"), momentum, innerTrackTimeReso);
+ histos.fill(HIST("h2dInnerTimeResTotalPionVsP"), momentum, innerTotalTimeReso);
+ }
+ if (ii == 3 && std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[3])->PdgCode()) {
+ histos.fill(HIST("h2dInnerTimeResTrackKaonVsP"), momentum, innerTrackTimeReso);
+ histos.fill(HIST("h2dInnerTimeResTotalKaonVsP"), momentum, innerTotalTimeReso);
+ }
+ if (ii == 4 && std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[4])->PdgCode()) {
+ histos.fill(HIST("h2dInnerTimeResTrackProtVsP"), momentum, innerTrackTimeReso);
+ histos.fill(HIST("h2dInnerTimeResTotalProtVsP"), momentum, innerTotalTimeReso);
+ }
+ }
+ if (doQAplots && trackLengthRecoOuterTOF > 0) {
+ float momentum = recoTrack.getP();
+ float pseudorapidity = recoTrack.getEta();
+ float transverse_momentum = momentum / std::cosh(pseudorapidity);
+ if (ii == 0 && std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[0])->PdgCode()) {
+ histos.fill(HIST("h2dOuterTimeResTrackElecVsP"), momentum, outerTrackTimeReso);
+ histos.fill(HIST("h2dOuterTimeResTotalElecVsP"), momentum, outerTotalTimeReso);
+ }
+ if (ii == 1 && std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[1])->PdgCode()) {
+ histos.fill(HIST("h2dOuterTimeResTrackMuonVsP"), momentum, outerTrackTimeReso);
+ histos.fill(HIST("h2dOuterTimeResTotalMuonVsP"), momentum, outerTotalTimeReso);
+ }
+ if (ii == 2 && std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[2])->PdgCode()) {
+ histos.fill(HIST("h2dOuterTimeResTrackPionVsP"), momentum, outerTrackTimeReso);
+ histos.fill(HIST("h2dOuterTimeResTotalPionVsP"), momentum, outerTotalTimeReso);
+
+ histos.fill(HIST("h2dRelativePtResolution"), transverse_momentum, 100.0 * pt_resolution / transverse_momentum);
+ histos.fill(HIST("h2dRelativeEtaResolution"), pseudorapidity, 100.0 * eta_resolution / (std::fabs(pseudorapidity) + 1e-6));
+ }
+ if (ii == 3 && std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[3])->PdgCode()) {
+ histos.fill(HIST("h2dOuterTimeResTrackKaonVsP"), momentum, outerTrackTimeReso);
+ histos.fill(HIST("h2dOuterTimeResTotalKaonVsP"), momentum, outerTotalTimeReso);
+ }
+ if (ii == 4 && std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[4])->PdgCode()) {
+ histos.fill(HIST("h2dOuterTimeResTrackProtVsP"), momentum, outerTrackTimeReso);
+ histos.fill(HIST("h2dOuterTimeResTotalProtVsP"), momentum, outerTotalTimeReso);
+ }
+ }
+ }
+
+ // Fixme: assumes dominant resolution effect is the TOF resolution
+ // and not the tracking itself. It's *probably* a fair assumption
+ // but it should be tested further! --> FIXED IN THIS VERSION
+ if (trackLengthInnerTOF > 0 && trackLengthRecoInnerTOF > 0)
+ nSigmaInnerTOF[ii] = deltaTimeInnerTOF[ii] / innerTotalTimeReso;
+ if (trackLengthOuterTOF > 0 && trackLengthRecoOuterTOF > 0)
+ nSigmaOuterTOF[ii] = deltaTimeOuterTOF[ii] / outerTotalTimeReso;
+ }
+
+ if (doQAplots) {
+ float momentum = recoTrack.getP();
+ if (trackLengthRecoInnerTOF > 0) {
+ if (std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[0])->PdgCode()) {
+ histos.fill(HIST("h2dInnerNsigmaTrueElecVsElecHypothesis"), momentum, nSigmaInnerTOF[0]);
+ histos.fill(HIST("h2dInnerNsigmaTrueElecVsMuonHypothesis"), momentum, nSigmaInnerTOF[1]);
+ histos.fill(HIST("h2dInnerNsigmaTrueElecVsPionHypothesis"), momentum, nSigmaInnerTOF[2]);
+ histos.fill(HIST("h2dInnerNsigmaTrueElecVsKaonHypothesis"), momentum, nSigmaInnerTOF[3]);
+ histos.fill(HIST("h2dInnerNsigmaTrueElecVsProtHypothesis"), momentum, nSigmaInnerTOF[4]);
+ }
+ if (std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[1])->PdgCode()) {
+ histos.fill(HIST("h2dInnerNsigmaTrueMuonVsElecHypothesis"), momentum, nSigmaInnerTOF[0]);
+ histos.fill(HIST("h2dInnerNsigmaTrueMuonVsMuonHypothesis"), momentum, nSigmaInnerTOF[1]);
+ histos.fill(HIST("h2dInnerNsigmaTrueMuonVsPionHypothesis"), momentum, nSigmaInnerTOF[2]);
+ histos.fill(HIST("h2dInnerNsigmaTrueMuonVsKaonHypothesis"), momentum, nSigmaInnerTOF[3]);
+ histos.fill(HIST("h2dInnerNsigmaTrueMuonVsProtHypothesis"), momentum, nSigmaInnerTOF[4]);
+ }
+ if (std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[2])->PdgCode()) {
+ histos.fill(HIST("h2dInnerNsigmaTruePionVsElecHypothesis"), momentum, nSigmaInnerTOF[0]);
+ histos.fill(HIST("h2dInnerNsigmaTruePionVsMuonHypothesis"), momentum, nSigmaInnerTOF[1]);
+ histos.fill(HIST("h2dInnerNsigmaTruePionVsPionHypothesis"), momentum, nSigmaInnerTOF[2]);
+ histos.fill(HIST("h2dInnerNsigmaTruePionVsKaonHypothesis"), momentum, nSigmaInnerTOF[3]);
+ histos.fill(HIST("h2dInnerNsigmaTruePionVsProtHypothesis"), momentum, nSigmaInnerTOF[4]);
+ }
+ if (std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[3])->PdgCode()) {
+ histos.fill(HIST("h2dInnerNsigmaTrueKaonVsElecHypothesis"), momentum, nSigmaInnerTOF[0]);
+ histos.fill(HIST("h2dInnerNsigmaTrueKaonVsMuonHypothesis"), momentum, nSigmaInnerTOF[1]);
+ histos.fill(HIST("h2dInnerNsigmaTrueKaonVsPionHypothesis"), momentum, nSigmaInnerTOF[2]);
+ histos.fill(HIST("h2dInnerNsigmaTrueKaonVsKaonHypothesis"), momentum, nSigmaInnerTOF[3]);
+ histos.fill(HIST("h2dInnerNsigmaTrueKaonVsProtHypothesis"), momentum, nSigmaInnerTOF[4]);
+ }
+ if (std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[4])->PdgCode()) {
+ histos.fill(HIST("h2dInnerNsigmaTrueProtVsElecHypothesis"), momentum, nSigmaInnerTOF[0]);
+ histos.fill(HIST("h2dInnerNsigmaTrueProtVsMuonHypothesis"), momentum, nSigmaInnerTOF[1]);
+ histos.fill(HIST("h2dInnerNsigmaTrueProtVsPionHypothesis"), momentum, nSigmaInnerTOF[2]);
+ histos.fill(HIST("h2dInnerNsigmaTrueProtVsKaonHypothesis"), momentum, nSigmaInnerTOF[3]);
+ histos.fill(HIST("h2dInnerNsigmaTrueProtVsProtHypothesis"), momentum, nSigmaInnerTOF[4]);
+ }
+ }
+ if (trackLengthRecoOuterTOF > 0) {
+ if (std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[0])->PdgCode()) {
+ histos.fill(HIST("h2dOuterNsigmaTrueElecVsElecHypothesis"), momentum, nSigmaOuterTOF[0]);
+ histos.fill(HIST("h2dOuterNsigmaTrueElecVsMuonHypothesis"), momentum, nSigmaOuterTOF[1]);
+ histos.fill(HIST("h2dOuterNsigmaTrueElecVsPionHypothesis"), momentum, nSigmaOuterTOF[2]);
+ histos.fill(HIST("h2dOuterNsigmaTrueElecVsKaonHypothesis"), momentum, nSigmaOuterTOF[3]);
+ histos.fill(HIST("h2dOuterNsigmaTrueElecVsProtHypothesis"), momentum, nSigmaOuterTOF[4]);
+ }
+ if (std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[1])->PdgCode()) {
+ histos.fill(HIST("h2dOuterNsigmaTrueMuonVsElecHypothesis"), momentum, nSigmaOuterTOF[0]);
+ histos.fill(HIST("h2dOuterNsigmaTrueMuonVsMuonHypothesis"), momentum, nSigmaOuterTOF[1]);
+ histos.fill(HIST("h2dOuterNsigmaTrueMuonVsPionHypothesis"), momentum, nSigmaOuterTOF[2]);
+ histos.fill(HIST("h2dOuterNsigmaTrueMuonVsKaonHypothesis"), momentum, nSigmaOuterTOF[3]);
+ histos.fill(HIST("h2dOuterNsigmaTrueMuonVsProtHypothesis"), momentum, nSigmaOuterTOF[4]);
+ }
+ if (std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[2])->PdgCode()) {
+ histos.fill(HIST("h2dOuterNsigmaTruePionVsElecHypothesis"), momentum, nSigmaOuterTOF[0]);
+ histos.fill(HIST("h2dOuterNsigmaTruePionVsMuonHypothesis"), momentum, nSigmaOuterTOF[1]);
+ histos.fill(HIST("h2dOuterNsigmaTruePionVsPionHypothesis"), momentum, nSigmaOuterTOF[2]);
+ histos.fill(HIST("h2dOuterNsigmaTruePionVsKaonHypothesis"), momentum, nSigmaOuterTOF[3]);
+ histos.fill(HIST("h2dOuterNsigmaTruePionVsProtHypothesis"), momentum, nSigmaOuterTOF[4]);
+ }
+ if (std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[3])->PdgCode()) {
+ histos.fill(HIST("h2dOuterNsigmaTrueKaonVsElecHypothesis"), momentum, nSigmaOuterTOF[0]);
+ histos.fill(HIST("h2dOuterNsigmaTrueKaonVsMuonHypothesis"), momentum, nSigmaOuterTOF[1]);
+ histos.fill(HIST("h2dOuterNsigmaTrueKaonVsPionHypothesis"), momentum, nSigmaOuterTOF[2]);
+ histos.fill(HIST("h2dOuterNsigmaTrueKaonVsKaonHypothesis"), momentum, nSigmaOuterTOF[3]);
+ histos.fill(HIST("h2dOuterNsigmaTrueKaonVsProtHypothesis"), momentum, nSigmaOuterTOF[4]);
+ }
+ if (std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(lpdg_array[4])->PdgCode()) {
+ histos.fill(HIST("h2dOuterNsigmaTrueProtVsElecHypothesis"), momentum, nSigmaOuterTOF[0]);
+ histos.fill(HIST("h2dOuterNsigmaTrueProtVsMuonHypothesis"), momentum, nSigmaOuterTOF[1]);
+ histos.fill(HIST("h2dOuterNsigmaTrueProtVsPionHypothesis"), momentum, nSigmaOuterTOF[2]);
+ histos.fill(HIST("h2dOuterNsigmaTrueProtVsKaonHypothesis"), momentum, nSigmaOuterTOF[3]);
+ histos.fill(HIST("h2dOuterNsigmaTrueProtVsProtHypothesis"), momentum, nSigmaOuterTOF[4]);
+ }
+ }
+ }
+
+ float deltaTrackLengthInnerTOF = abs(trackLengthInnerTOF - trackLengthRecoInnerTOF);
+ if (trackLengthInnerTOF > 0 && trackLengthRecoInnerTOF > 0) {
+ histos.fill(HIST("h2dDeltaTrackLengthInnerVsPt"), o2track.getPt(), deltaTrackLengthInnerTOF);
+ }
+ float deltaTrackLengthOuterTOF = abs(trackLengthOuterTOF - trackLengthRecoOuterTOF);
+ if (trackLengthOuterTOF > 0 && trackLengthRecoOuterTOF > 0) {
+ histos.fill(HIST("h2dDeltaTrackLengthOuterVsPt"), o2track.getPt(), deltaTrackLengthInnerTOF);
+ }
+
+ // Sigmas have been fully calculated. Please populate the NSigma helper table (once per track)
+ upgradeTof(nSigmaInnerTOF[0], nSigmaInnerTOF[1], nSigmaInnerTOF[2], nSigmaInnerTOF[3], nSigmaInnerTOF[4], trackLengthInnerTOF, trackLengthRecoInnerTOF, deltaTrackLengthInnerTOF,
+ nSigmaOuterTOF[0], nSigmaOuterTOF[1], nSigmaOuterTOF[2], nSigmaOuterTOF[3], nSigmaOuterTOF[4], trackLengthOuterTOF, trackLengthRecoOuterTOF, deltaTrackLengthOuterTOF);
+ }
+ }
+};
+
+WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
+{
+ return WorkflowSpec{adaptAnalysisTask(cfgc)};
+}
diff --git a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx
new file mode 100644
index 00000000000..531d1d088cb
--- /dev/null
+++ b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx
@@ -0,0 +1,552 @@
+// 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 onTheFlyTracker.cxx
+///
+/// \brief LUT-based on-the-fly analysis task-level tracking
+///
+/// This task allows for the calculation of aod::collisions and aod::Tracks in a synthetic manner,
+/// smearing MC particles with very configurable settings. This will allow for the usage of
+/// custom LUTs (obtained through separate studies) and the subsequent estimate of the performance
+/// of a future detector even in very statistics-hungry analyses.
+///
+/// \author David Dobrigkeit Chinellato, UNICAMP
+/// \author Nicolò Jacazio , UniBo
+/// \author Roberto Preghenella preghenella@bo.infn.it
+///
+
+#include
+
+#include
+
+#include "Framework/AnalysisDataModel.h"
+#include "Framework/AnalysisTask.h"
+#include "Framework/runDataProcessing.h"
+#include "Framework/HistogramRegistry.h"
+#include
+#include "Framework/O2DatabasePDGPlugin.h"
+#include "Common/DataModel/TrackSelectionTables.h"
+#include "ReconstructionDataFormats/DCA.h"
+#include "DetectorsBase/Propagator.h"
+#include "DataFormatsParameters/GRPMagField.h"
+#include "DetectorsVertexing/PVertexer.h"
+#include "DetectorsVertexing/PVertexerHelpers.h"
+#include "SimulationDataFormat/InteractionSampler.h"
+#include "Field/MagneticField.h"
+
+#include "ALICE3/Core/DelphesO2TrackSmearer.h"
+#include "ALICE3/DataModel/collisionAlice3.h"
+#include "ALICE3/DataModel/tracksAlice3.h"
+
+using namespace o2;
+using namespace o2::framework;
+
+struct OnTheFlyTracker {
+ Produces collisions;
+ Produces collLabels;
+ Produces tracksPar;
+ Produces tracksParExtension;
+ Produces tracksParCov;
+ Produces tracksParCovExtension;
+ Produces tracksLabels;
+ Produces tracksDCA;
+ Produces collisionsAlice3;
+ Produces TracksAlice3;
+
+ // optionally produced, empty (to be tuned later)
+ Produces tracksExtra; // base table, extend later
+ Produces trackSelection;
+ Produces trackSelectionExtension;
+
+ Configurable magneticField{"magneticField", 20.0f, "magnetic field in kG"};
+ Configurable maxEta{"maxEta", 1.5, "maximum eta to consider viable"};
+ Configurable multEtaRange{"multEtaRange", 0.8, "eta range to compute the multiplicity"};
+ Configurable minPt{"minPt", 0.1, "minimum pt to consider viable"};
+ Configurable enableLUT{"enableLUT", false, "Enable track smearing"};
+ Configurable enableNucleiSmearing{"enableNucleiSmearing", false, "Enable smearing of nuclei"};
+ Configurable enablePrimaryVertexing{"enablePrimaryVertexing", true, "Enable primary vertexing"};
+ Configurable interpolateLutEfficiencyVsNch{"interpolateLutEfficiencyVsNch", true, "interpolate LUT efficiency as f(Nch)"};
+
+ Configurable populateTracksDCA{"populateTracksDCA", true, "populate TracksDCA table"};
+ Configurable populateTracksExtra{"populateTracksExtra", false, "populate TracksExtra table (legacy)"};
+ Configurable populateTrackSelection{"populateTrackSelection", false, "populate TrackSelection table (legacy)"};
+
+ Configurable processUnreconstructedTracks{"processUnreconstructedTracks", false, "process (smear) unreco-ed tracks"};
+ Configurable doExtraQA{"doExtraQA", false, "do extra 2D QA plots"};
+ Configurable extraQAwithoutDecayDaughters{"extraQAwithoutDecayDaughters", false, "remove decay daughters from qa plots (yes/no)"};
+
+ Configurable lutEl{"lutEl", "lutCovm.el.dat", "LUT for electrons"};
+ Configurable lutMu{"lutMu", "lutCovm.mu.dat", "LUT for muons"};
+ Configurable lutPi{"lutPi", "lutCovm.pi.dat", "LUT for pions"};
+ Configurable lutKa{"lutKa", "lutCovm.ka.dat", "LUT for kaons"};
+ Configurable lutPr{"lutPr", "lutCovm.pr.dat", "LUT for protons"};
+ Configurable lutDe{"lutDe", "lutCovm.de.dat", "LUT for deuterons"};
+ Configurable lutTr{"lutTr", "lutCovm.tr.dat", "LUT for tritons"};
+ Configurable lutHe3{"lutHe3", "lutCovm.he3.dat", "LUT for Helium-3"};
+
+ ConfigurableAxis axisMomentum{"axisMomentum", {VARIABLE_WIDTH, 0.0f, 0.1f, 0.2f, 0.3f, 0.4f, 0.5f, 0.6f, 0.7f, 0.8f, 0.9f, 1.0f, 1.1f, 1.2f, 1.3f, 1.4f, 1.5f, 1.6f, 1.7f, 1.8f, 1.9f, 2.0f, 2.2f, 2.4f, 2.6f, 2.8f, 3.0f, 3.2f, 3.4f, 3.6f, 3.8f, 4.0f, 4.4f, 4.8f, 5.2f, 5.6f, 6.0f, 6.5f, 7.0f, 7.5f, 8.0f, 9.0f, 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 17.0f, 19.0f, 21.0f, 23.0f, 25.0f, 30.0f, 35.0f, 40.0f, 50.0f}, "#it{p} (GeV/#it{c})"};
+ ConfigurableAxis axisNVertices{"axisNVertices", {20, -0.5, 19.5}, "N_{vertices}"};
+ ConfigurableAxis axisMultiplicity{"axisMultiplicity", {100, -0.5, 99.5}, "N_{contributors}"};
+ ConfigurableAxis axisVertexZ{"axisVertexZ", {40, -20, 20}, "vertex Z (cm)"};
+ ConfigurableAxis axisDCA{"axisDCA", {400, -200, 200}, "DCA (#mum)"};
+ ConfigurableAxis axisX{"axisX", {250, -50, 200}, "track X (cm)"};
+
+ using PVertex = o2::dataformats::PrimaryVertex;
+
+ // Class to hold the track information for the O2 vertexing
+ class TrackAlice3 : public o2::track::TrackParCov
+ {
+ using TimeEst = o2::dataformats::TimeStampWithError;
+
+ public:
+ TrackAlice3() = default;
+ ~TrackAlice3() = default;
+ TrackAlice3(const TrackAlice3& src) = default;
+ TrackAlice3(const o2::track::TrackParCov& src, const int64_t label, const float t = 0, const float te = 1, bool decayDauInput = false) : o2::track::TrackParCov(src), mcLabel{label}, timeEst{t, te}, isDecayDau(decayDauInput) {}
+ const TimeEst& getTimeMUS() const { return timeEst; }
+ int64_t mcLabel;
+ TimeEst timeEst; ///< time estimate in ns
+ bool isDecayDau;
+ };
+
+ // necessary for particle charges
+ Service pdgDB;
+
+ // for handling basic QA histograms if requested
+ HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject};
+
+ o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrNONE;
+
+ // Track smearer
+ o2::delphes::DelphesO2TrackSmearer mSmearer;
+
+ // For processing and vertexing
+ std::vector tracksAlice3;
+ std::vector ghostTracksAlice3;
+ std::vector bcData;
+ o2::steer::InteractionSampler irSampler;
+ o2::vertexing::PVertexer vertexer;
+
+ void init(o2::framework::InitContext& initContext)
+ {
+ if (enableLUT) {
+ std::map mapPdgLut;
+ const char* lutElChar = lutEl->c_str();
+ const char* lutMuChar = lutMu->c_str();
+ const char* lutPiChar = lutPi->c_str();
+ const char* lutKaChar = lutKa->c_str();
+ const char* lutPrChar = lutPr->c_str();
+
+ LOGF(info, "Will load electron lut file ..: %s", lutElChar);
+ LOGF(info, "Will load muon lut file ......: %s", lutMuChar);
+ LOGF(info, "Will load pion lut file ......: %s", lutPiChar);
+ LOGF(info, "Will load kaon lut file ......: %s", lutKaChar);
+ LOGF(info, "Will load proton lut file ....: %s", lutPrChar);
+
+ mapPdgLut.insert(std::make_pair(11, lutElChar));
+ mapPdgLut.insert(std::make_pair(13, lutMuChar));
+ mapPdgLut.insert(std::make_pair(211, lutPiChar));
+ mapPdgLut.insert(std::make_pair(321, lutKaChar));
+ mapPdgLut.insert(std::make_pair(2212, lutPrChar));
+
+ if (enableNucleiSmearing) {
+ const char* lutDeChar = lutDe->c_str();
+ const char* lutTrChar = lutTr->c_str();
+ const char* lutHe3Char = lutHe3->c_str();
+ mapPdgLut.insert(std::make_pair(1000010020, lutDeChar));
+ mapPdgLut.insert(std::make_pair(1000010030, lutTrChar));
+ mapPdgLut.insert(std::make_pair(1000020030, lutHe3Char));
+ }
+ for (auto e : mapPdgLut) {
+ if (!mSmearer.loadTable(e.first, e.second)) {
+ LOG(fatal) << "Having issue with loading the LUT " << e.first << " " << e.second;
+ }
+ }
+ // interpolate efficiencies if requested to do so
+ mSmearer.interpolateEfficiency(static_cast(interpolateLutEfficiencyVsNch));
+
+ // smear un-reco'ed tracks if asked to do so
+ mSmearer.skipUnreconstructed(static_cast(!processUnreconstructedTracks));
+ }
+
+ // Basic QA
+ histos.add("hPtGenerated", "hPtGenerated", kTH1F, {axisMomentum});
+ histos.add("hPtGeneratedEl", "hPtGeneratedEl", kTH1F, {axisMomentum});
+ histos.add("hPtGeneratedPi", "hPtGeneratedPi", kTH1F, {axisMomentum});
+ histos.add("hPtGeneratedKa", "hPtGeneratedKa", kTH1F, {axisMomentum});
+ histos.add("hPtGeneratedPr", "hPtGeneratedPr", kTH1F, {axisMomentum});
+ histos.add("hPtReconstructed", "hPtReconstructed", kTH1F, {axisMomentum});
+ histos.add("hPtReconstructedEl", "hPtReconstructedEl", kTH1F, {axisMomentum});
+ histos.add("hPtReconstructedPi", "hPtReconstructedPi", kTH1F, {axisMomentum});
+ histos.add("hPtReconstructedKa", "hPtReconstructedKa", kTH1F, {axisMomentum});
+ histos.add("hPtReconstructedPr", "hPtReconstructedPr", kTH1F, {axisMomentum});
+
+ // Collision QA
+ histos.add("hPVz", "hPVz", kTH1F, {axisVertexZ});
+ histos.add("hLUTMultiplicity", "hLUTMultiplicity", kTH1F, {axisMultiplicity});
+ histos.add("hSimMultiplicity", "hSimMultiplicity", kTH1F, {axisMultiplicity});
+ histos.add("hRecoMultiplicity", "hRecoMultiplicity", kTH1F, {axisMultiplicity});
+
+ if (doExtraQA) {
+ histos.add("h2dVerticesVsContributors", "h2dVerticesVsContributors", kTH2F, {axisMultiplicity, axisNVertices});
+ histos.add("hRecoVsSimMultiplicity", "hRecoVsSimMultiplicity", kTH2F, {axisMultiplicity, axisMultiplicity});
+ histos.add("h2dDCAxy", "h2dDCAxy", kTH2F, {axisMomentum, axisDCA});
+
+ histos.add("hSimTrackX", "hSimTrackX", kTH1F, {axisX});
+ histos.add("hRecoTrackX", "hRecoTrackX", kTH1F, {axisX});
+ histos.add("hTrackXatDCA", "hTrackXatDCA", kTH1F, {axisX});
+ }
+
+ LOGF(info, "Initializing magnetic field to value: %.3f kG", static_cast(magneticField));
+ o2::parameters::GRPMagField grpmag;
+ grpmag.setFieldUniformity(true);
+ grpmag.setL3Current(30000.f * (magneticField / 5.0f));
+ auto field = grpmag.getNominalL3Field();
+ o2::base::Propagator::initFieldFromGRP(&grpmag);
+
+ auto fieldInstance = static_cast(TGeoGlobalMagField::Instance()->GetField());
+ if (!fieldInstance) {
+ LOGF(fatal, "Failed to set up magnetic field! Stopping now!");
+ }
+
+ // Cross-check
+ LOGF(info, "Check field at (0, 0, 0): %.1f kG, nominal: %.1f", static_cast(fieldInstance->GetBz(0, 0, 0)), static_cast(field));
+
+ LOGF(info, "Initializing empty material cylinder LUT - could be better in the future");
+ o2::base::MatLayerCylSet* lut = new o2::base::MatLayerCylSet();
+ lut->addLayer(200, 200.1, 2, 1.0f, 100.0f);
+ LOGF(info, "MatLayerCylSet::optimizePhiSlices()");
+ lut->optimizePhiSlices();
+ LOGF(info, "Setting lut now...");
+ o2::base::Propagator::Instance()->setMatLUT(lut);
+
+ irSampler.setInteractionRate(100);
+ irSampler.setFirstIR(o2::InteractionRecord(0, 0));
+ irSampler.init();
+
+ vertexer.setValidateWithIR(kFALSE);
+ vertexer.setBunchFilling(irSampler.getBunchFilling());
+ vertexer.init();
+ }
+
+ /// Function to convert a McParticle into a perfect Track
+ /// \param particle the particle to convert (mcParticle)
+ /// \param o2track the address of the resulting TrackParCov
+ template
+ void convertMCParticleToO2Track(McParticleType& particle, o2::track::TrackParCov& o2track)
+ {
+ auto pdgInfo = pdgDB->GetParticle(particle.pdgCode());
+ int charge = 0;
+ if (pdgInfo != nullptr) {
+ charge = pdgInfo->Charge() / 3;
+ }
+ std::array params;
+ std::array covm = {0.};
+ float s, c, x;
+ o2::math_utils::sincos(particle.phi(), s, c);
+ o2::math_utils::rotateZInv(particle.vx(), particle.vy(), x, params[0], s, c);
+ params[1] = particle.vz();
+ params[2] = 0.; // since alpha = phi
+ auto theta = 2. * std::atan(std::exp(-particle.eta()));
+ params[3] = 1. / std::tan(theta);
+ params[4] = charge / particle.pt();
+
+ // Initialize TrackParCov in-place
+ new (&o2track)(o2::track::TrackParCov)(x, particle.phi(), params, covm);
+ }
+
+ float dNdEta = 0.f; // Charged particle multiplicity to use in the efficiency evaluation
+ void process(aod::McCollision const& mcCollision, aod::McParticles const& mcParticles)
+ {
+ tracksAlice3.clear();
+ ghostTracksAlice3.clear();
+ bcData.clear();
+
+ o2::dataformats::DCA dcaInfo;
+ o2::dataformats::VertexBase vtx;
+
+ // generate collision time
+ auto ir = irSampler.generateCollisionTime();
+
+ // First we compute the number of charged particles in the event
+ dNdEta = 0.f;
+ for (const auto& mcParticle : mcParticles) {
+ if (std::abs(mcParticle.eta()) > multEtaRange) {
+ continue;
+ }
+ if (!mcParticle.isPhysicalPrimary()) {
+ continue;
+ }
+ const auto pdg = std::abs(mcParticle.pdgCode());
+ if (pdg != kElectron && pdg != kMuonMinus && pdg != kPiPlus && pdg != kKPlus && pdg != kProton) {
+ continue;
+ }
+ const auto& pdgInfo = pdgDB->GetParticle(mcParticle.pdgCode());
+ if (!pdgInfo) {
+ LOG(warning) << "PDG code " << mcParticle.pdgCode() << " not found in the database";
+ continue;
+ }
+ if (pdgInfo->Charge() == 0) {
+ continue;
+ }
+ dNdEta += 1.f;
+ }
+
+ dNdEta /= (multEtaRange * 2.0f);
+ uint32_t multiplicityCounter = 0;
+ histos.fill(HIST("hLUTMultiplicity"), dNdEta);
+
+ for (const auto& mcParticle : mcParticles) {
+ if (!mcParticle.isPhysicalPrimary()) {
+ continue;
+ }
+ const auto pdg = std::abs(mcParticle.pdgCode());
+ if (pdg != kElectron && pdg != kMuonMinus && pdg != kPiPlus && pdg != kKPlus && pdg != kProton) {
+ continue;
+ }
+ if (std::fabs(mcParticle.eta()) > maxEta) {
+ continue;
+ }
+
+ histos.fill(HIST("hPtGenerated"), mcParticle.pt());
+ if (TMath::Abs(mcParticle.pdgCode()) == 11)
+ histos.fill(HIST("hPtGeneratedEl"), mcParticle.pt());
+ if (TMath::Abs(mcParticle.pdgCode()) == 211)
+ histos.fill(HIST("hPtGeneratedPi"), mcParticle.pt());
+ if (TMath::Abs(mcParticle.pdgCode()) == 321)
+ histos.fill(HIST("hPtGeneratedKa"), mcParticle.pt());
+ if (TMath::Abs(mcParticle.pdgCode()) == 2212)
+ histos.fill(HIST("hPtGeneratedPr"), mcParticle.pt());
+
+ if (mcParticle.pt() < minPt) {
+ continue;
+ }
+
+ bool isDecayDaughter = false;
+ if (mcParticle.getProcess() == 4)
+ isDecayDaughter = true;
+
+ multiplicityCounter++;
+ o2::track::TrackParCov trackParCov;
+ convertMCParticleToO2Track(mcParticle, trackParCov);
+
+ if (doExtraQA) {
+ histos.fill(HIST("hSimTrackX"), trackParCov.getX());
+ }
+
+ bool reconstructed = mSmearer.smearTrack(trackParCov, mcParticle.pdgCode(), dNdEta);
+ if (!reconstructed && !processUnreconstructedTracks) {
+ continue;
+ }
+ if (TMath::IsNaN(trackParCov.getZ())) {
+ // capture rare smearing mistakes / corrupted tracks
+ continue;
+ }
+
+ // Base QA (note: reco pT here)
+ histos.fill(HIST("hPtReconstructed"), trackParCov.getPt());
+ if (TMath::Abs(mcParticle.pdgCode()) == 11)
+ histos.fill(HIST("hPtReconstructedEl"), mcParticle.pt());
+ if (TMath::Abs(mcParticle.pdgCode()) == 211)
+ histos.fill(HIST("hPtReconstructedPi"), mcParticle.pt());
+ if (TMath::Abs(mcParticle.pdgCode()) == 321)
+ histos.fill(HIST("hPtReconstructedKa"), mcParticle.pt());
+ if (TMath::Abs(mcParticle.pdgCode()) == 2212)
+ histos.fill(HIST("hPtReconstructedPr"), mcParticle.pt());
+
+ if (doExtraQA) {
+ histos.fill(HIST("hRecoTrackX"), trackParCov.getX());
+ }
+
+ // populate vector with track if we reco-ed it
+ const float t = (ir.timeInBCNS + gRandom->Gaus(0., 100.)) * 1e-3;
+ if (reconstructed) {
+ tracksAlice3.push_back(TrackAlice3{trackParCov, mcParticle.globalIndex(), t, 100.f * 1e-3, isDecayDaughter});
+ } else {
+ ghostTracksAlice3.push_back(TrackAlice3{trackParCov, mcParticle.globalIndex(), t, 100.f * 1e-3, isDecayDaughter});
+ }
+ }
+
+ // *+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*
+ // Calculate primary vertex with tracks from this collision
+ // data preparation
+ o2::vertexing::PVertex primaryVertex;
+
+ if (enablePrimaryVertexing) {
+ std::vector lblTracks;
+ std::vector vertices;
+ std::vector vertexTrackIDs;
+ std::vector v2tRefs;
+ std::vector lblVtx;
+ lblVtx.emplace_back(mcCollision.globalIndex(), 1);
+ std::vector idxVec; // store IDs
+
+ idxVec.reserve(tracksAlice3.size());
+ for (unsigned i = 0; i < tracksAlice3.size(); i++) {
+ lblTracks.emplace_back(tracksAlice3[i].mcLabel, mcCollision.globalIndex(), 1, false);
+ idxVec.emplace_back(i, o2::dataformats::GlobalTrackID::ITS); // let's say ITS
+ }
+
+ // Calculate vertices
+ const int n_vertices = vertexer.process(tracksAlice3, // track array
+ idxVec,
+ gsl::span{bcData},
+ vertices,
+ vertexTrackIDs,
+ v2tRefs,
+ gsl::span{lblTracks},
+ lblVtx);
+
+ if (n_vertices < 1) {
+ return; // primary vertex not reconstructed
+ }
+
+ // Find largest vertex
+ int largestVertex = 0;
+ for (Int_t iv = 1; iv < n_vertices; iv++) {
+ if (vertices[iv].getNContributors() > vertices[largestVertex].getNContributors()) {
+ largestVertex = iv;
+ }
+ }
+ primaryVertex = vertices[largestVertex];
+ if (doExtraQA) {
+ histos.fill(HIST("h2dVerticesVsContributors"), primaryVertex.getNContributors(), n_vertices);
+ }
+ } else {
+ primaryVertex.setXYZ(mcCollision.posX(), mcCollision.posY(), mcCollision.posZ());
+ }
+ // *+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*
+
+ // debug / informational
+ histos.fill(HIST("hSimMultiplicity"), multiplicityCounter);
+ histos.fill(HIST("hRecoMultiplicity"), tracksAlice3.size());
+ histos.fill(HIST("hPVz"), primaryVertex.getZ());
+
+ if (doExtraQA) {
+ histos.fill(HIST("hRecoVsSimMultiplicity"), multiplicityCounter, tracksAlice3.size());
+ }
+
+ // *+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*
+ // populate collisions
+ collisions(-1, // BC is irrelevant in synthetic MC tests for now, could be adjusted in future
+ primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ(),
+ primaryVertex.getSigmaX2(), primaryVertex.getSigmaXY(), primaryVertex.getSigmaY2(),
+ primaryVertex.getSigmaXZ(), primaryVertex.getSigmaYZ(), primaryVertex.getSigmaZ2(),
+ 0, primaryVertex.getChi2(), primaryVertex.getNContributors(),
+ 0, 0);
+ collLabels(mcCollision.globalIndex(), 0);
+ collisionsAlice3(dNdEta);
+ // *+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*
+
+ // *+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*
+ // populate tracks
+ for (const auto& trackParCov : tracksAlice3) {
+ // Fixme: collision index could be changeable
+ aod::track::TrackTypeEnum trackType = aod::track::Track;
+
+ if (populateTracksDCA) {
+ float dcaXY = 1e+10, dcaZ = 1e+10;
+ o2::track::TrackParCov trackParametrization(trackParCov);
+ if (trackParametrization.propagateToDCA(primaryVertex, magneticField, &dcaInfo)) {
+ dcaXY = dcaInfo.getY();
+ dcaZ = dcaInfo.getZ();
+ }
+ if (doExtraQA && (!extraQAwithoutDecayDaughters || (extraQAwithoutDecayDaughters && !trackParCov.isDecayDau))) {
+ histos.fill(HIST("h2dDCAxy"), trackParametrization.getPt(), dcaXY * 1e+4); // in microns, please
+ histos.fill(HIST("hTrackXatDCA"), trackParametrization.getX());
+ }
+ tracksDCA(dcaXY, dcaZ);
+ }
+
+ tracksPar(collisions.lastIndex(), trackType, trackParCov.getX(), trackParCov.getAlpha(), trackParCov.getY(), trackParCov.getZ(), trackParCov.getSnp(), trackParCov.getTgl(), trackParCov.getQ2Pt());
+ tracksParExtension(trackParCov.getPt(), trackParCov.getP(), trackParCov.getEta(), trackParCov.getPhi());
+
+ // TODO do we keep the rho as 0? Also the sigma's are duplicated information
+ tracksParCov(std::sqrt(trackParCov.getSigmaY2()), std::sqrt(trackParCov.getSigmaZ2()), std::sqrt(trackParCov.getSigmaSnp2()),
+ std::sqrt(trackParCov.getSigmaTgl2()), std::sqrt(trackParCov.getSigma1Pt2()), 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
+ tracksParCovExtension(trackParCov.getSigmaY2(), trackParCov.getSigmaZY(), trackParCov.getSigmaZ2(), trackParCov.getSigmaSnpY(),
+ trackParCov.getSigmaSnpZ(), trackParCov.getSigmaSnp2(), trackParCov.getSigmaTglY(), trackParCov.getSigmaTglZ(), trackParCov.getSigmaTglSnp(),
+ trackParCov.getSigmaTgl2(), trackParCov.getSigma1PtY(), trackParCov.getSigma1PtZ(), trackParCov.getSigma1PtSnp(), trackParCov.getSigma1PtTgl(),
+ trackParCov.getSigma1Pt2());
+ tracksLabels(trackParCov.mcLabel, 0);
+
+ // populate extra tables if required to do so
+ if (populateTracksExtra) {
+ tracksExtra(0.0f, (uint32_t)0, (uint8_t)0, (uint8_t)0,
+ (int8_t)0, (int8_t)0, (uint8_t)0, (uint8_t)0,
+ 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f,
+ 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f);
+ }
+ if (populateTrackSelection) {
+ trackSelection((uint8_t)0, false, false, false, false, false, false);
+ trackSelectionExtension(false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false);
+ }
+ TracksAlice3(true);
+ }
+ // populate ghost tracks
+ for (const auto& trackParCov : ghostTracksAlice3) {
+ // Fixme: collision index could be changeable
+ aod::track::TrackTypeEnum trackType = aod::track::Track;
+
+ if (populateTracksDCA) {
+ float dcaXY = 1e+10, dcaZ = 1e+10;
+ o2::track::TrackParCov trackParametrization(trackParCov);
+ if (trackParametrization.propagateToDCA(primaryVertex, magneticField, &dcaInfo)) {
+ dcaXY = dcaInfo.getY();
+ dcaZ = dcaInfo.getZ();
+ }
+ if (doExtraQA && (!extraQAwithoutDecayDaughters || (extraQAwithoutDecayDaughters && !trackParCov.isDecayDau))) {
+ histos.fill(HIST("h2dDCAxy"), trackParametrization.getPt(), dcaXY * 1e+4); // in microns, please
+ histos.fill(HIST("hTrackXatDCA"), trackParametrization.getX());
+ }
+ tracksDCA(dcaXY, dcaZ);
+ }
+
+ tracksPar(collisions.lastIndex(), trackType, trackParCov.getX(), trackParCov.getAlpha(), trackParCov.getY(), trackParCov.getZ(), trackParCov.getSnp(), trackParCov.getTgl(), trackParCov.getQ2Pt());
+ tracksParExtension(trackParCov.getPt(), trackParCov.getP(), trackParCov.getEta(), trackParCov.getPhi());
+
+ // TODO do we keep the rho as 0? Also the sigma's are duplicated information
+ tracksParCov(std::sqrt(trackParCov.getSigmaY2()), std::sqrt(trackParCov.getSigmaZ2()), std::sqrt(trackParCov.getSigmaSnp2()),
+ std::sqrt(trackParCov.getSigmaTgl2()), std::sqrt(trackParCov.getSigma1Pt2()), 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
+ tracksParCovExtension(trackParCov.getSigmaY2(), trackParCov.getSigmaZY(), trackParCov.getSigmaZ2(), trackParCov.getSigmaSnpY(),
+ trackParCov.getSigmaSnpZ(), trackParCov.getSigmaSnp2(), trackParCov.getSigmaTglY(), trackParCov.getSigmaTglZ(), trackParCov.getSigmaTglSnp(),
+ trackParCov.getSigmaTgl2(), trackParCov.getSigma1PtY(), trackParCov.getSigma1PtZ(), trackParCov.getSigma1PtSnp(), trackParCov.getSigma1PtTgl(),
+ trackParCov.getSigma1Pt2());
+ tracksLabels(trackParCov.mcLabel, 0);
+
+ // populate extra tables if required to do so
+ if (populateTracksExtra) {
+ tracksExtra(0.0f, (uint32_t)0, (uint8_t)0, (uint8_t)0,
+ (int8_t)0, (int8_t)0, (uint8_t)0, (uint8_t)0,
+ 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f,
+ 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f);
+ }
+ if (populateTrackSelection) {
+ trackSelection((uint8_t)0, false, false, false, false, false, false);
+ trackSelectionExtension(false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false);
+ }
+ TracksAlice3(false);
+ }
+ }
+};
+
+/// Extends TracksExtra if necessary
+struct onTheFlyTrackerInitializer {
+ Spawns tracksExtra;
+ void init(InitContext const&) {}
+};
+
+WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
+{
+ return WorkflowSpec{
+ adaptAnalysisTask(cfgc),
+ adaptAnalysisTask(cfgc)};
+}
diff --git a/ALICE3/TableProducer/alice3-decayfinder.cxx b/ALICE3/TableProducer/alice3-decayfinder.cxx
new file mode 100644
index 00000000000..04f8284e3c4
--- /dev/null
+++ b/ALICE3/TableProducer/alice3-decayfinder.cxx
@@ -0,0 +1,598 @@
+// 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.
+//
+// *+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*
+// Decay finder task for ALICE 3
+// *+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*
+//
+// Uses specific ALICE 3 PID and performance for studying
+// HF decays. Work in progress: use at your own risk!
+//
+
+#include
+#include
+#include
+#include