diff --git a/analysis/detectors/sts/CMakeLists.txt b/analysis/detectors/sts/CMakeLists.txt
index 6b4a0ac3ee0192e1da1e1aa94c4cc2aaac055a34..ecb7c91bb259de480f64f7729a8b260a607c43df 100644
--- a/analysis/detectors/sts/CMakeLists.txt
+++ b/analysis/detectors/sts/CMakeLists.txt
@@ -7,12 +7,29 @@ set(INCLUDE_DIRECTORIES
set(SRCS
CbmStsWkn.cxx
+ CbmCutId.cxx
+ CbmCutMap.cxx
+ CbmStsUtils.cxx
+ CbmStsAnaBase.cxx
+ CbmStsAnalysis.cxx
+ CbmSpillCheck.cxx
+ CbmStsChannelQA.cxx
+ CbmStsTimeCal.cxx
+ CbmStsHitAna.cxx
+ CbmStsCorrelation.cxx
+ CbmStsRecoBeamSpot.cxx
+ CbmDcaVertexFinder.cxx
+ CbmStsEfficiency.cxx
+ CbmStsResolution.cxx
+ CbmEventVertexDca.cxx
)
-
set(LIBRARY_NAME CbmStsAna)
set(LINKDEF ${LIBRARY_NAME}LinkDef.h)
set(PUBLIC_DEPENDENCIES
+ Algo
+ CbmSimBase
+ CbmStsBase
ROOT::Core
)
@@ -22,4 +39,11 @@ set(PRIVATE_DEPENDENCIES
ROOT::MathCore
)
+install(FILES ${CMAKE_CURRENT_SOURCE_DIR}/CbmCut.h
+ DESTINATION include)
+
generate_cbm_library()
+
+If(GTEST_FOUND)
+ add_subdirectory(tests)
+EndIf()
diff --git a/analysis/detectors/sts/CbmCut.h b/analysis/detectors/sts/CbmCut.h
new file mode 100644
index 0000000000000000000000000000000000000000..fac886230e776742e86bce9ffaed7f4f1f5e89f6
--- /dev/null
+++ b/analysis/detectors/sts/CbmCut.h
@@ -0,0 +1,90 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+#ifndef CBMCUT_H
+#define CBMCUT_H
+
+#include "CbmCutId.h"
+
+#include <type_traits>
+
+template<typename T>
+class CbmCut {
+ public:
+ /**
+ * @brief Default constructor
+ * Any check done to a non-configured object CbmCut will pass
+ */
+ CbmCut() = default;
+
+ /**
+ * @brief Add a CbmTarget object to the list of targets with a key as trg_name.
+ * @param min Minimum value for the check
+ * @param max Maximum value for the check
+ * When this constructor is created, check will be done for each boundary
+ */
+ CbmCut(T min, T max)
+ {
+ SetMin(min);
+ SetMax(max);
+ }
+
+ /**
+ * @brief Check if the T is inside the configured ranges [min, max].
+ * @param a Object for which operators <= >= == is implemented
+ * @return true or false depending on two possible general cases
+ * Case 1
+ * -----------------########## TRUE ##########-----------------
+ * @MIN @MAX
+ *
+ * Case 2
+ * #### TRUE ####--------------------------------#### TRUE ####
+ * @MAX @MIN
+ */
+ bool Check(T a) const
+ {
+ bool check_min_status = fMinState ? a >= min_ : true;
+ bool check_max_status = fMaxState ? a <= max_ : true;
+ if (fMinState && fMaxState && (max_ < min_)) return check_min_status || check_max_status;
+ return check_min_status && check_max_status;
+ }
+
+ void SetMin(T val)
+ {
+ min_ = val;
+ fMinState = true;
+ }
+ void SetMax(T val)
+ {
+ max_ = val;
+ fMaxState = true;
+ }
+ void SetRange(T min_val, T max_val)
+ {
+ SetMin(min_val);
+ SetMax(max_val);
+ }
+
+ friend std::ostream& operator<<(std::ostream& out, const CbmCut<T>& obj)
+ {
+ std::string min_msg = obj.fMinState ? std::to_string(obj.min_) : "none";
+ std::string max_msg = obj.fMaxState ? std::to_string(obj.max_) : "none";
+ out << "[Min: " << min_msg << " , Max: " << max_msg << "]";
+ return out;
+ }
+
+ T GetMin() const { return min_; }
+ T GetMax() const { return max_; }
+ bool GetMinState() const { return fMinState; }
+ bool GetMaxState() const { return fMaxState; }
+
+
+ private:
+ bool fMinState{false};
+ bool fMaxState{false};
+ T min_;
+ T max_;
+};
+
+#endif
diff --git a/analysis/detectors/sts/CbmCutId.cxx b/analysis/detectors/sts/CbmCutId.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..e39fb3eacb1a6b125fcd12ffae0e7705f80dfecf
--- /dev/null
+++ b/analysis/detectors/sts/CbmCutId.cxx
@@ -0,0 +1,20 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+#include "CbmCutId.h"
+
+std::string ToString(CbmCutId id)
+{
+ auto result = cut_id_to_str_map.find(id);
+ if (result == cut_id_to_str_map.end()) {
+ return "NotExist";
+ }
+ return result->second;
+};
+
+std::ostream& operator<<(std::ostream& strm, const CbmCutId& cut_id)
+{
+ strm << ToString(cut_id);
+ return strm;
+}
diff --git a/analysis/detectors/sts/CbmCutId.h b/analysis/detectors/sts/CbmCutId.h
new file mode 100644
index 0000000000000000000000000000000000000000..747b982f8ab931fc7b122c9303cef232341a8147
--- /dev/null
+++ b/analysis/detectors/sts/CbmCutId.h
@@ -0,0 +1,257 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+#ifndef CBMCUTID_H
+#define CBMCUTID_H
+
+#include <iostream>
+#include <unordered_map>
+
+/** \enum CbmCutId
+ * \brief Enumeration of cut identifiers for various observables.
+ */
+enum class CbmCutId : ushort
+{
+ // Digi
+ kBmonDigiChannel,
+ kBmonDigiSide,
+
+ kMvdDigiChannel,
+ kMvdDigiCharge,
+ kMvdDigiTime,
+
+ kStsDigiChannel,
+ kStsDigiCharge,
+ kStsDigiTime,
+
+ kRichDigiChannel,
+ kRichDigiCharge,
+ kRichDigiTime,
+
+ kMuchigiChannel,
+ kMuchigiCharge,
+ kMuchigiTime,
+
+ kTrdDigiChannel,
+ kTrdDigiCharge,
+ kTrdDigiTime,
+
+ kTofDigiChannel,
+ kTofDigiCharge,
+ kTofDigiTime,
+
+ // Cluster
+ kStsClusterAddress,
+ kStsClusterTime,
+ kStsClusterCharge,
+ kStsClusterSize,
+ kStsClusterPosition,
+
+ // Hit
+ kMvdHitAddress,
+ kMvdHitTime,
+ kMvdHitCharge,
+ kMvdHitX,
+ kMvdHitY,
+ kMvdHitZ,
+
+ kStsHitAddress,
+ kStsHitTime,
+ kStsHitCharge,
+ kStsHitQasym,
+ kStsHitX,
+ kStsHitY,
+ kStsHitZ,
+
+ kRichHitAddress,
+ kRichHitTime,
+ kRichHitCharge,
+ kRichHitX,
+ kRichHitY,
+ kRichHitZ,
+
+ kMuchHitAddress,
+ kMuchHitTime,
+ kMuchHitCharge,
+ kMuchHitX,
+ kMuchHitY,
+ kMuchHitZ,
+
+ kTrdHitAddress,
+ kTrdHitTime,
+ kTrdHitCharge,
+ kTrdHitX,
+ kTrdHitY,
+ kTrdHitZ,
+
+ kTofHitAddress,
+ kTofHitTime,
+ kTofHitCharge,
+ kTofHitX,
+ kTofHitY,
+ kTofHitZ,
+
+ // Event
+ kEventNofMvdHit,
+ kEventNofStsHit,
+ kEventNofRichHit,
+ kEventNofRichRing,
+ kEventNofMuchPixelHit,
+ kEventNofMuchStrawHit,
+ kEventNofTrdHit,
+ kEventNofTofHit,
+ kEventNofGlobalTrack,
+
+ // Detector tracks
+ kMvdTrackChi2,
+ kMvdTrackSize,
+
+ kStsTrackChi2,
+ kStsTrackSize,
+
+ kRichTrackChi2,
+ kRichTrackSize,
+
+ kMuchTrackChi2,
+ kMuchTrackSize,
+
+ kTrdTrackChi2,
+ kTrdTrackSize,
+
+ kTofTrackChi2,
+ kTofTrackSize,
+
+ // Global tracks
+ kGlobalTrackChi2,
+ kGlobalTrackPval,
+ kGlobalTrackSize,
+ kGlobalTrackLength,
+
+ kGlobalTrackMvdChi2,
+ kGlobalTrackMvdSize,
+
+ kGlobalTrackStsChi2,
+ kGlobalTrackStsSize,
+
+ kGlobalTrackRichChi2,
+ kGlobalTrackRichSize,
+
+ kGlobalTrackMuchChi2,
+ kGlobalTrackMuchSize,
+
+ kGlobalTrackTrdChi2,
+ kGlobalTrackTrdSize,
+
+ kGlobalTrackTofChi2,
+ kGlobalTrackTofSize
+};
+
+static const std::unordered_map<CbmCutId, std::string> cut_id_to_str_map = {
+ {CbmCutId::kBmonDigiChannel, "kBmonDigiChannel"},
+ {CbmCutId::kBmonDigiSide, "kBmonDigiSide"},
+ {CbmCutId::kMvdDigiChannel, "kMvdDigiChannel"},
+ {CbmCutId::kMvdDigiCharge, "kMvdDigiCharge"},
+ {CbmCutId::kMvdDigiTime, "kMvdDigiTime"},
+ {CbmCutId::kStsDigiChannel, "kStsDigiChannel"},
+ {CbmCutId::kStsDigiCharge, "kStsDigiCharge"},
+ {CbmCutId::kStsDigiTime, "kStsDigiTime"},
+ {CbmCutId::kRichDigiChannel, "kRichDigiChannel"},
+ {CbmCutId::kRichDigiCharge, "kRichDigiCharge"},
+ {CbmCutId::kRichDigiTime, "kRichDigiTime"},
+ {CbmCutId::kMuchigiChannel, "kMuchigiChannel"},
+ {CbmCutId::kMuchigiCharge, "kMuchigiCharge"},
+ {CbmCutId::kMuchigiTime, "kMuchigiTime"},
+ {CbmCutId::kTrdDigiChannel, "kTrdDigiChannel"},
+ {CbmCutId::kTrdDigiCharge, "kTrdDigiCharge"},
+ {CbmCutId::kTrdDigiTime, "kTrdDigiTime"},
+ {CbmCutId::kTofDigiChannel, "kTofDigiChannel"},
+ {CbmCutId::kTofDigiCharge, "kTofDigiCharge"},
+ {CbmCutId::kTofDigiTime, "kTofDigiTime"},
+ {CbmCutId::kStsClusterAddress, "kStsClusterAddress"},
+ {CbmCutId::kStsClusterTime, "kStsClusterTime"},
+ {CbmCutId::kStsClusterCharge, "kStsClusterCharge"},
+ {CbmCutId::kStsClusterSize, "kStsClusterSize"},
+ {CbmCutId::kStsClusterPosition, "kStsClusterPosition"},
+ {CbmCutId::kMvdHitAddress, "kMvdHitAddress"},
+ {CbmCutId::kMvdHitTime, "kMvdHitTime"},
+ {CbmCutId::kMvdHitCharge, "kMvdHitCharge"},
+ {CbmCutId::kMvdHitX, "kMvdHitX"},
+ {CbmCutId::kMvdHitY, "kMvdHitY"},
+ {CbmCutId::kMvdHitZ, "kMvdHitZ"},
+ {CbmCutId::kStsHitAddress, "kStsHitAddress"},
+ {CbmCutId::kStsHitTime, "kStsHitTime"},
+ {CbmCutId::kStsHitCharge, "kStsHitCharge"},
+ {CbmCutId::kStsHitQasym, "kStsHitQasym"},
+ {CbmCutId::kStsHitX, "kStsHitX"},
+ {CbmCutId::kStsHitY, "kStsHitY"},
+ {CbmCutId::kStsHitZ, "kStsHitZ"},
+ {CbmCutId::kRichHitAddress, "kRichHitAddress"},
+ {CbmCutId::kRichHitTime, "kRichHitTime"},
+ {CbmCutId::kRichHitCharge, "kRichHitCharge"},
+ {CbmCutId::kRichHitX, "kRichHitX"},
+ {CbmCutId::kRichHitY, "kRichHitY"},
+ {CbmCutId::kRichHitZ, "kRichHitZ"},
+ {CbmCutId::kMuchHitAddress, "kMuchHitAddress"},
+ {CbmCutId::kMuchHitTime, "kMuchHitTime"},
+ {CbmCutId::kMuchHitCharge, "kMuchHitCharge"},
+ {CbmCutId::kMuchHitX, "kMuchHitX"},
+ {CbmCutId::kMuchHitY, "kMuchHitY"},
+ {CbmCutId::kMuchHitZ, "kMuchHitZ"},
+ {CbmCutId::kTrdHitAddress, "kTrdHitAddress"},
+ {CbmCutId::kTrdHitTime, "kTrdHitTime"},
+ {CbmCutId::kTrdHitCharge, "kTrdHitCharge"},
+ {CbmCutId::kTrdHitX, "kTrdHitX"},
+ {CbmCutId::kTrdHitY, "kTrdHitY"},
+ {CbmCutId::kTrdHitZ, "kTrdHitZ"},
+ {CbmCutId::kTofHitAddress, "kTofHitAddress"},
+ {CbmCutId::kTofHitTime, "kTofHitTime"},
+ {CbmCutId::kTofHitCharge, "kTofHitCharge"},
+ {CbmCutId::kTofHitX, "kTofHitX"},
+ {CbmCutId::kTofHitY, "kTofHitY"},
+ {CbmCutId::kTofHitZ, "kTofHitZ"},
+ {CbmCutId::kEventNofMvdHit, "kEventNofMvdHit"},
+ {CbmCutId::kEventNofStsHit, "kEventNofStsHit"},
+ {CbmCutId::kEventNofRichHit, "kEventNofRichHit"},
+ {CbmCutId::kEventNofRichRing, "kEventNofRichRing"},
+ {CbmCutId::kEventNofMuchPixelHit, "kEventNofMuchPixelHit"},
+ {CbmCutId::kEventNofMuchStrawHit, "kEventNofMuchStrawHit"},
+ {CbmCutId::kEventNofTrdHit, "kEventNofTrdHit"},
+ {CbmCutId::kEventNofTofHit, "kEventNofTofHit"},
+ {CbmCutId::kEventNofGlobalTrack, "kEventNofGlobalTrack"},
+ {CbmCutId::kMvdTrackChi2, "kMvdTrackChi2"},
+ {CbmCutId::kMvdTrackSize, "kMvdTrackSize"},
+ {CbmCutId::kStsTrackChi2, "kStsTrackChi2"},
+ {CbmCutId::kStsTrackSize, "kStsTrackSize"},
+ {CbmCutId::kRichTrackChi2, "kRichTrackChi2"},
+ {CbmCutId::kRichTrackSize, "kRichTrackSize"},
+ {CbmCutId::kMuchTrackChi2, "kMuchTrackChi2"},
+ {CbmCutId::kMuchTrackSize, "kMuchTrackSize"},
+ {CbmCutId::kTrdTrackChi2, "kTrdTrackChi2"},
+ {CbmCutId::kTrdTrackSize, "kTrdTrackSize"},
+ {CbmCutId::kTofTrackChi2, "kTofTrackChi2"},
+ {CbmCutId::kTofTrackSize, "kTofTrackSize"},
+ {CbmCutId::kGlobalTrackChi2, "kGlobalTrackChi2"},
+ {CbmCutId::kGlobalTrackPval, "kGlobalTrackPval"},
+ {CbmCutId::kGlobalTrackSize, "kGlobalTrackSize"},
+ {CbmCutId::kGlobalTrackLength, "kGlobalTrackLength"},
+ {CbmCutId::kGlobalTrackMvdChi2, "kGlobalTrackMvdChi2"},
+ {CbmCutId::kGlobalTrackMvdSize, "kGlobalTrackMvdSize"},
+ {CbmCutId::kGlobalTrackStsChi2, "kGlobalTrackStsChi2"},
+ {CbmCutId::kGlobalTrackStsSize, "kGlobalTrackStsSize"},
+ {CbmCutId::kGlobalTrackRichChi2, "kGlobalTrackRichChi2"},
+ {CbmCutId::kGlobalTrackRichSize, "kGlobalTrackRichSize"},
+ {CbmCutId::kGlobalTrackMuchChi2, "kGlobalTrackMuchChi2"},
+ {CbmCutId::kGlobalTrackMuchSize, "kGlobalTrackMuchSize"},
+ {CbmCutId::kGlobalTrackTrdChi2, "kGlobalTrackTrdChi2"},
+ {CbmCutId::kGlobalTrackTrdSize, "kGlobalTrackTrdSize"},
+ {CbmCutId::kGlobalTrackTofChi2, "kGlobalTrackTofChi2"},
+ {CbmCutId::kGlobalTrackTofSize, "kGlobalTrackTofSize"}};
+
+/** \brief Convert CbmCutId to a string representation.
+ *
+ * \param cutId The CbmCutId to convert.
+ * \return String representation of the CbmCutId.
+ */
+std::string ToString(CbmCutId);
+#endif
diff --git a/analysis/detectors/sts/CbmCutMap.cxx b/analysis/detectors/sts/CbmCutMap.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..b055814bee61f166301554f7c2c881f02f551eca
--- /dev/null
+++ b/analysis/detectors/sts/CbmCutMap.cxx
@@ -0,0 +1,49 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+#include "CbmCutMap.h"
+
+CbmCut<float>* CbmCutMap::AddCbmCut(CbmCutId id) { return &fCbmCuts[id]; }
+
+bool CbmCutMap::Check(CbmCutId id, double value)
+{
+ if (!fCbmCuts.count(id)) {
+ fFailedCounter[id] = 0;
+ return true;
+ }
+ bool check_result = fCbmCuts[id].Check(value);
+ fFailedCounter[id] = check_result ? fFailedCounter[id] + 1 : fFailedCounter[id];
+ return check_result;
+}
+
+bool CbmCutMap::CheckStsHit(CbmStsHit* hit, TClonesArray* clu_array = nullptr)
+{
+ if (hit == nullptr) return false;
+
+ if (clu_array != nullptr) {
+ if (!Check(CbmCutId::kStsHitCharge, cbm_sts_utils::GetHitCharge(hit, clu_array))) {
+ return false;
+ }
+ if (!Check(CbmCutId::kStsHitQasym, cbm_sts_utils::GetHitChargeAsy(hit, clu_array))) {
+ return false;
+ }
+ }
+ return Check(CbmCutId::kStsHitTime, hit->GetTime()) && Check(CbmCutId::kStsHitX, hit->GetX())
+ && Check(CbmCutId::kStsHitY, hit->GetY()) && Check(CbmCutId::kStsHitZ, hit->GetZ());
+}
+
+bool CbmCutMap::CheckEvent(CbmEvent* evt)
+{
+ if (evt == nullptr) return false;
+
+ return Check(CbmCutId::kEventNofMvdHit, evt->GetNofData(ECbmDataType::kMvdHit))
+ && Check(CbmCutId::kEventNofStsHit, evt->GetNofData(ECbmDataType::kStsHit))
+ && Check(CbmCutId::kEventNofRichHit, evt->GetNofData(ECbmDataType::kRichHit))
+ && Check(CbmCutId::kEventNofRichRing, evt->GetNofData(ECbmDataType::kRichRing))
+ && Check(CbmCutId::kEventNofMuchPixelHit, evt->GetNofData(ECbmDataType::kMuchPixelHit))
+ && Check(CbmCutId::kEventNofMuchStrawHit, evt->GetNofData(ECbmDataType::kMuchStrawHit))
+ && Check(CbmCutId::kEventNofTrdHit, evt->GetNofData(ECbmDataType::kTrdHit))
+ && Check(CbmCutId::kEventNofTofHit, evt->GetNofData(ECbmDataType::kTofHit))
+ && Check(CbmCutId::kEventNofGlobalTrack, evt->GetNofData(ECbmDataType::kGlobalTrack));
+}
diff --git a/analysis/detectors/sts/CbmCutMap.h b/analysis/detectors/sts/CbmCutMap.h
new file mode 100644
index 0000000000000000000000000000000000000000..513e9b9e98bf73a8c77470176cb3fbb566571542
--- /dev/null
+++ b/analysis/detectors/sts/CbmCutMap.h
@@ -0,0 +1,89 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+#ifndef CBMCUTMAP_H
+#define CBMCUTMAP_H
+
+#include "CbmCut.h"
+#include "CbmEvent.h"
+#include "CbmStsHit.h"
+#include "CbmStsUtils.h"
+
+#include <TClonesArray.h>
+
+/**
+ * @class CbmCutMap
+ * @brief Manages a map of cuts associated with their IDs.
+ */
+class CbmCutMap {
+ public:
+ CbmCutMap() = default;
+ ~CbmCutMap() = default;
+
+ /**
+ * @brief Get the map of cuts.
+ * @return The map of cuts.
+ */
+ std::unordered_map<CbmCutId, CbmCut<float>> GetMap() const { return fCbmCuts; }
+
+ /**
+ * @brief Add a new cut to the map.
+ * @param id The ID of the cut.
+ * @return Pointer to the added CbmCut.
+ */
+ CbmCut<float>* AddCbmCut(CbmCutId id);
+
+ /**
+ * @brief Check if a value passes the cut with the given ID.
+ * @param id The ID of the cut to check.
+ * @param value The value to check against the cut.
+ * @return True if the value passes the cut, false otherwise.
+ */
+ bool Check(CbmCutId id, double value);
+
+ /**
+ * @brief Check if a CbmStsHit passes the cuts.
+ * @param hit Pointer to the CbmStsHit to check.
+ * @param array Pointer to the TClonesArray containing CbmStsCluster (for charge retrieving).
+ * @return True if the hit passes the cuts, false otherwise.
+ */
+ bool CheckStsHit(CbmStsHit*, TClonesArray*);
+
+ /**
+ * @brief Check if a CbmEvent passes the cuts.
+ * @param evt Pointer to the CbmEvent to check.
+ * @return True if the event passes the cuts, false otherwise.
+ */
+ bool CheckEvent(CbmEvent* evt);
+
+ /**
+ * @brief Overloaded stream insertion operator for CbmCutMap.
+ * @param out Output stream.
+ * @param obj CbmCutMap object to output.
+ * @return Reference to the output stream.
+ */
+ friend std::ostream& operator<<(std::ostream& out, const CbmCutMap& obj)
+ {
+ for (auto& [id, cut] : obj.GetMap()) {
+ out << ToString(id) << "\t" << cut << std::endl;
+ }
+ return out;
+ }
+
+ /**
+ * @brief Print the cuts and failed pass counters.
+ */
+ void Print()
+ {
+ for (auto& [id, cut] : fCbmCuts) {
+ std::cout << int(id) << ": " << cut << "\tFailed: " << fFailedCounter[id] << std::endl;
+ }
+ }
+
+ private:
+ std::unordered_map<CbmCutId, CbmCut<float>> fCbmCuts;
+ std::unordered_map<CbmCutId, unsigned long int> fFailedCounter;
+};
+
+#endif
diff --git a/analysis/detectors/sts/CbmDcaVertexFinder.cxx b/analysis/detectors/sts/CbmDcaVertexFinder.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..bd394f0466e00ca2f162adb82b0dbb4d8ca5ce81
--- /dev/null
+++ b/analysis/detectors/sts/CbmDcaVertexFinder.cxx
@@ -0,0 +1,127 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+#include "CbmDcaVertexFinder.h"
+
+#include "CbmGlobalTrack.h"
+#include "CbmVertex.h"
+
+#include <Logger.h>
+
+CbmDcaVertexFinder::CbmDcaVertexFinder()
+ : fMaxDca(0.1)
+ , fInputTracks(std::vector<CbmGlobalTrack*>())
+ , fNbPairs(0)
+ , fQA(std::nullopt)
+{
+ LOG(debug) << "Creating CbmDcaVertexFinder ...";
+ fCovMatrix.ResizeTo(3, 3);
+}
+
+CbmDcaVertexFinder::CbmDcaVertexFinder(double max_dca)
+ : fMaxDca(max_dca)
+ , fInputTracks(std::vector<CbmGlobalTrack*>())
+ , fNbPairs(0)
+ , fQA(std::nullopt)
+{
+ LOG(debug) << "Creating CbmDcaVertexFinder ...";
+ fCovMatrix.ResizeTo(3, 3);
+}
+
+/** \brief Default constructor
+ * \param tracks
+*/
+CbmDcaVertexFinder::CbmDcaVertexFinder(const std::vector<CbmGlobalTrack*> tracks)
+ : fMaxDca(0.1)
+ , fInputTracks(tracks)
+ , fNbPairs(0)
+ , fQA(std::nullopt)
+{
+ LOG(debug) << "Creating CbmDcaVertexFinder ...";
+ fCovMatrix.ResizeTo(3, 3);
+}
+
+/** \brief Default constructor
+ * \param tracks
+*/
+CbmDcaVertexFinder::CbmDcaVertexFinder(const std::vector<CbmGlobalTrack*> tracks, const double max_dca)
+ : fMaxDca(max_dca)
+ , fInputTracks(tracks)
+ , fNbPairs(0)
+ , fQA(std::nullopt)
+{
+ LOG(debug) << "Creating CbmDcaVertexFinder ...";
+ fCovMatrix.ResizeTo(3, 3);
+}
+
+void CbmDcaVertexFinder::SetTracks(const std::vector<CbmGlobalTrack*> tracks)
+{
+ fInputTracks.clear();
+ fInputTracks = tracks;
+}
+
+/** \brief find vertex using input tracks
+ */
+std::optional<CbmVertex> CbmDcaVertexFinder::FindVertex()
+{
+ // Reset the number of track pair used
+ fNbPairs = 0;
+ int n_of_trk = fInputTracks.size();
+ LOG(debug) << "- PCA - Find event vertex using CbmGlobalTracks: " << n_of_trk;
+ TVector3 vtx;
+ for (int trk_i_idx = 0; trk_i_idx < n_of_trk - 1; trk_i_idx++) {
+ for (int trk_j_idx = trk_i_idx + 1; trk_j_idx < n_of_trk; trk_j_idx++) {
+ auto pca = FindPca(fInputTracks[trk_i_idx], fInputTracks[trk_j_idx]);
+ if (pca.has_value() && pca->d_trk < fMaxDca) {
+ TVector3 pca_i_j = pca->point;
+ vtx += pca_i_j;
+ fNbPairs++;
+
+ if (fQA.has_value()) {
+ fQA->pca_y_vs_x->Fill(pca_i_j[0], pca_i_j[1]);
+ fQA->pca_x_vs_z->Fill(pca_i_j[2], pca_i_j[0]);
+ fQA->pca_y_vs_z->Fill(pca_i_j[2], pca_i_j[1]);
+ }
+ }
+ }
+ }
+ if (!fNbPairs) return std::nullopt;
+
+ vtx *= 1. / fNbPairs;
+
+ // WARNING LINE
+ return CbmVertex("EventVertex", "EventVertex", vtx[0], vtx[1], vtx[2], 0, 0, fNbPairs, fCovMatrix);
+}
+
+std::optional<CbmDcaVertexFinder::PCA> CbmDcaVertexFinder::FindPca(CbmGlobalTrack* trk_i, CbmGlobalTrack* trk_j)
+{
+ TVector3 r1(trk_i->GetParamFirst()->GetX(), trk_i->GetParamFirst()->GetY(), trk_i->GetParamFirst()->GetZ());
+ TVector3 e1(trk_i->GetParamFirst()->GetTx(), trk_i->GetParamFirst()->GetTy(), 1);
+
+ TVector3 r2(trk_j->GetParamFirst()->GetX(), trk_j->GetParamFirst()->GetY(), trk_j->GetParamFirst()->GetZ());
+ TVector3 e2(trk_j->GetParamFirst()->GetTx(), trk_j->GetParamFirst()->GetTy(), 1);
+
+ TVector3 n = e1.Cross(e2); // Directional vector for segment of the pca
+
+ float nn = n.Dot(n);
+ if (nn != 0) { // Non-Parallel tracks
+
+ TVector3 e1n = e1.Cross(n);
+ TVector3 e2n = e1.Cross(n);
+ TVector3 r21 = r2 - r1;
+
+ float t1 = e2n.Dot(r21) / nn;
+ float t2 = e1n.Dot(r21) / nn;
+
+ TVector3 p1 = r1 + t1 * e1; // Closest point on track 1
+ TVector3 p2 = r2 + t2 * e2; // Closest point on track 2
+ TVector3 p21 = p2 - p1;
+
+ TVector3 point = 0.5 * (p1 + p2);
+ double d_trk = 0.5 * p21.Mag();
+
+ return std::make_optional<PCA>({point, d_trk});
+ }
+ return std::nullopt;
+}
diff --git a/analysis/detectors/sts/CbmDcaVertexFinder.h b/analysis/detectors/sts/CbmDcaVertexFinder.h
new file mode 100644
index 0000000000000000000000000000000000000000..f3830f17e315bfd58e2d98ddd05cb64fdb7326fe
--- /dev/null
+++ b/analysis/detectors/sts/CbmDcaVertexFinder.h
@@ -0,0 +1,75 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+#ifndef CBMDCAVERTEXFINDER_H
+#define CBMDCAVERTEXFINDER_H
+
+#include "TH2.h"
+#include "TMatrixFSym.h"
+#include "TVector3.h"
+
+#include <optional>
+class CbmVertex;
+class CbmGlobalTrack;
+
+/** \class CbmDcaVertexFinder
+ * \brief Estimates a common vertex from multiple straight GlobalTracks
+ *
+ * The algorithm consist in averaging the PCA of all valid pair of tracks.
+ * A valid pair is defined but a cut in the maximum value for DCA
+ */
+class CbmDcaVertexFinder {
+ public:
+ struct Qa {
+ std::shared_ptr<TH2> pca_y_vs_x;
+ std::shared_ptr<TH2> pca_x_vs_z;
+ std::shared_ptr<TH2> pca_y_vs_z;
+ };
+
+ /** Holds PCA analysis for a pair of straight tracks */
+ struct PCA {
+ TVector3 point;
+ double d_trk;
+ };
+
+ CbmDcaVertexFinder();
+
+ /** \brief Constructor setting maximum DCA for track pair */
+ CbmDcaVertexFinder(double max_dca);
+
+ /** \brief Constructor using input set of track */
+ CbmDcaVertexFinder(const std::vector<CbmGlobalTrack*> tracks);
+
+ /** \brief Constructor using input set of track, and maximum DCA for track pair */
+ CbmDcaVertexFinder(const std::vector<CbmGlobalTrack*> tracks, const double max_dca);
+
+ /** \brief Set input track */
+ void SetTracks(const std::vector<CbmGlobalTrack*> tracks);
+
+ ~CbmDcaVertexFinder() = default;
+
+ /** \brief Calculate the Point of Closest Approach of two straight tracks
+ * if tracks are parallel it return a nullopt
+ * \param trk_i, trk_j pointer to CbmGlobal tracks
+ */
+ std::optional<CbmDcaVertexFinder::PCA> FindPca(CbmGlobalTrack* trk_i, CbmGlobalTrack* trk_j);
+
+ /** \brief Run algorithm to find vertex
+ * \return std::optional<CbmVertex>:
+ * - std::nullopt if algorithm fails or no vertex was found under present configuration
+ */
+ std::optional<CbmVertex> FindVertex();
+
+ void EnableQa(Qa qa) { fQA = std::make_optional<Qa>(qa); }
+
+ private:
+ const double fMaxDca;
+ std::vector<CbmGlobalTrack*> fInputTracks;
+ double fNbPairs;
+ std::optional<Qa> fQA;
+
+ TMatrixFSym fCovMatrix;
+};
+
+#endif
diff --git a/analysis/detectors/sts/CbmEventVertexDca.cxx b/analysis/detectors/sts/CbmEventVertexDca.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..f5cc0bb53aadda5ec371a3a40ee74c55b548c09a
--- /dev/null
+++ b/analysis/detectors/sts/CbmEventVertexDca.cxx
@@ -0,0 +1,354 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+#include "CbmEventVertexDca.h"
+
+#include "CbmStsAddress.h"
+
+
+void CbmEventVertexDca::BookHistograms()
+{
+ std::string h_name;
+
+ // ---- Binning configuration ----
+ // 3D Vertex
+ double vertex_dx = 0.01;
+ double vertex_dy = 0.01;
+ double vertex_dz = 0.01;
+ double vertex_x_min = -35 + 0.5 * vertex_dx;
+ double vertex_x_max = +35 + 0.5 * vertex_dx;
+ double vertex_y_min = -35 + 0.5 * vertex_dy;
+ double vertex_y_max = +35 + 0.5 * vertex_dy;
+ double vertex_z_min = -65 + 0.5 * vertex_dz;
+ double vertex_z_max = +65 + 0.5 * vertex_dz;
+
+ auto x_vtx_binning =
+ cbm_sts_utils::HBinning{uint32_t((vertex_x_max - vertex_x_min) / vertex_dx), vertex_x_min, vertex_x_max};
+ auto y_vtx_binning =
+ cbm_sts_utils::HBinning{uint32_t((vertex_y_max - vertex_y_min) / vertex_dy), vertex_y_min, vertex_y_max};
+ auto z_vtx_binning =
+ cbm_sts_utils::HBinning{uint32_t((vertex_z_max - vertex_z_min) / vertex_dz), vertex_z_min, vertex_z_max};
+
+ // DCA
+ double dca_dx = .001; // 10 um
+ double dca_min = -2.5 - 0.5 * dca_dx;
+ double dca_max = +2.5 + 0.5 * dca_dx;
+ auto r_dca_binning = cbm_sts_utils::HBinning{uint32_t((dca_max - dca_min) / dca_dx), dca_min, dca_max};
+ // ---- --------------------- ----
+
+ // ---- 3D Vertex 2D projections ----
+ h_name = "vertex_y_vs_x";
+ fH2D[h_name] =
+ std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), x_vtx_binning.n_of_bins, x_vtx_binning.x_min,
+ x_vtx_binning.x_max, y_vtx_binning.n_of_bins, y_vtx_binning.x_min, y_vtx_binning.x_max);
+ fH2D[h_name]->GetXaxis()->SetTitle("Vertex_{X} [cm]");
+ fH2D[h_name]->GetYaxis()->SetTitle("Vertex_{Y} [cm]");
+
+ h_name = "vertex_x_vs_z";
+ fH2D[h_name] =
+ std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), z_vtx_binning.n_of_bins, z_vtx_binning.x_min,
+ z_vtx_binning.x_max, x_vtx_binning.n_of_bins, x_vtx_binning.x_min, x_vtx_binning.x_max);
+ fH2D[h_name]->GetXaxis()->SetTitle("Vertex_{Z} [cm]");
+ fH2D[h_name]->GetYaxis()->SetTitle("Vertex_{X} [cm]");
+
+ h_name = "vertex_y_vs_z";
+ fH2D[h_name] =
+ std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), z_vtx_binning.n_of_bins, z_vtx_binning.x_min,
+ z_vtx_binning.x_max, y_vtx_binning.n_of_bins, y_vtx_binning.x_min, y_vtx_binning.x_max);
+ fH2D[h_name]->GetXaxis()->SetTitle("Vertex_{Z} [cm]");
+ fH2D[h_name]->GetYaxis()->SetTitle("Vertex_{Y} [cm]");
+ // ---- ------------------------ ----
+
+ // ---- PCA 2D projections ----
+ if (fVertexFinder) {
+ h_name = "pca_y_vs_x";
+ fH2DShared[h_name] =
+ std::make_shared<TH2D>(h_name.c_str(), h_name.c_str(), x_vtx_binning.n_of_bins, x_vtx_binning.x_min,
+ x_vtx_binning.x_max, y_vtx_binning.n_of_bins, y_vtx_binning.x_min, y_vtx_binning.x_max);
+ fH2DShared[h_name]->GetXaxis()->SetTitle("Vertex_{X} [cm]");
+ fH2DShared[h_name]->GetYaxis()->SetTitle("Vertex_{Y} [cm]");
+
+ h_name = "pca_x_vs_z";
+ fH2DShared[h_name] =
+ std::make_shared<TH2D>(h_name.c_str(), h_name.c_str(), z_vtx_binning.n_of_bins, z_vtx_binning.x_min,
+ z_vtx_binning.x_max, x_vtx_binning.n_of_bins, x_vtx_binning.x_min, x_vtx_binning.x_max);
+ fH2DShared[h_name]->GetXaxis()->SetTitle("Vertex_{Z} [cm]");
+ fH2DShared[h_name]->GetYaxis()->SetTitle("Vertex_{X} [cm]");
+
+ h_name = "pca_y_vs_z";
+ fH2DShared[h_name] =
+ std::make_shared<TH2D>(h_name.c_str(), h_name.c_str(), z_vtx_binning.n_of_bins, z_vtx_binning.x_min,
+ z_vtx_binning.x_max, y_vtx_binning.n_of_bins, y_vtx_binning.x_min, y_vtx_binning.x_max);
+ fH2DShared[h_name]->GetXaxis()->SetTitle("Vertex_{Z} [cm]");
+ fH2DShared[h_name]->GetYaxis()->SetTitle("Vertex_{Y} [cm]");
+
+ std::cout << " - DcaVertexFinder - Enabling QA" << std::endl;
+ fVertexFinder->EnableQa({fH2DShared["pca_y_vs_x"], fH2DShared["pca_x_vs_z"], fH2DShared["pca_y_vs_z"]});
+ }
+ // ---- ------------------------ ----
+
+ // ---- DCA ----
+ for (const char* di : {"x", "y", "z"}) {
+ h_name = Form("dca_d%s_vs_x", di);
+ fH2D[h_name] =
+ std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), r_dca_binning.n_of_bins, r_dca_binning.x_min,
+ r_dca_binning.x_max, x_vtx_binning.n_of_bins, x_vtx_binning.x_min, x_vtx_binning.x_max);
+ fH2D[h_name]->GetYaxis()->SetTitle("X [cm]");
+ fH2D[h_name]->GetXaxis()->SetTitle(Form("DCA_{%s} [cm]", di));
+
+ h_name = Form("dca_d%s_vs_y", di);
+ fH2D[h_name] =
+ std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), r_dca_binning.n_of_bins, r_dca_binning.x_min,
+ r_dca_binning.x_max, y_vtx_binning.n_of_bins, y_vtx_binning.x_min, y_vtx_binning.x_max);
+ fH2D[h_name]->GetYaxis()->SetTitle("Y [cm]");
+ fH2D[h_name]->GetXaxis()->SetTitle(Form("DCA_{%s} [cm]", di));
+
+ h_name = Form("dca_d%s_vs_z", di);
+ fH2D[h_name] =
+ std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), r_dca_binning.n_of_bins, r_dca_binning.x_min,
+ r_dca_binning.x_max, z_vtx_binning.n_of_bins, z_vtx_binning.x_min, z_vtx_binning.x_max);
+ fH2D[h_name]->GetYaxis()->SetTitle(Form("Z [cm]"));
+ fH2D[h_name]->GetXaxis()->SetTitle(Form("DCA_{%s} [cm]", di));
+ }
+ // ---- ------------------ ----
+
+ // ---- Residual with MC ----
+ if (input_has_mc) {
+ LOG(info) << " - INFO - Enabling MC comparison";
+ for (const char* di : {"x", "y", "z"}) {
+ h_name = Form("res_%s_vs_x", di);
+ fH2D[h_name] =
+ std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), r_dca_binning.n_of_bins, r_dca_binning.x_min,
+ r_dca_binning.x_max, x_vtx_binning.n_of_bins, x_vtx_binning.x_min, x_vtx_binning.x_max);
+ fH2D[h_name]->GetYaxis()->SetTitle("X [cm]");
+ fH2D[h_name]->GetXaxis()->SetTitle(Form("%s_{rec} - %s_{MC} [cm]", di, di));
+
+ h_name = Form("res_%s_vs_y", di);
+ fH2D[h_name] =
+ std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), r_dca_binning.n_of_bins, r_dca_binning.x_min,
+ r_dca_binning.x_max, y_vtx_binning.n_of_bins, y_vtx_binning.x_min, y_vtx_binning.x_max);
+ fH2D[h_name]->GetYaxis()->SetTitle("Y [cm]");
+ fH2D[h_name]->GetXaxis()->SetTitle(Form("%s_{rec} - %s_{MC} [cm]", di, di));
+
+ h_name = Form("res_%s_vs_z", di);
+ fH2D[h_name] =
+ std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), r_dca_binning.n_of_bins, r_dca_binning.x_min,
+ r_dca_binning.x_max, z_vtx_binning.n_of_bins, z_vtx_binning.x_min, z_vtx_binning.x_max);
+ fH2D[h_name]->GetYaxis()->SetTitle(Form("Z [cm]"));
+ fH2D[h_name]->GetXaxis()->SetTitle(Form("%s_{rec} - %s_{MC} [cm]", di, di));
+ }
+ }
+ // ---- ------------------------- ----
+
+ // Track multiplicity
+ h_name = "ca_track_multiplicity_all";
+ fH1D[h_name] = std::make_unique<TH1D>(h_name.c_str(), h_name.c_str(), 50, 0, 50);
+ fH1D[h_name]->GetXaxis()->SetTitle("N_{CATrack}");
+ fH1D[h_name]->GetYaxis()->SetTitle("Entries");
+
+ h_name = "ca_track_multiplicity_sel";
+ fH1D[h_name] = std::make_unique<TH1D>(h_name.c_str(), h_name.c_str(), 50, 0, 50);
+ fH1D[h_name]->GetXaxis()->SetTitle("N_{CATrack}");
+ fH1D[h_name]->GetYaxis()->SetTitle("Entries");
+}
+
+
+void CbmEventVertexDca::FinishTask() { SaveToFile(); }
+
+void CbmEventVertexDca::CheckVertex()
+{
+ fH1D["ca_track_multiplicity_sel"]->Fill(fGlbTrks.size());
+
+ fVertexFinder->SetTracks(fGlbTrks);
+ const std::optional<CbmVertex> vertex = fVertexFinder->FindVertex();
+
+ if (vertex.has_value()) {
+
+ TVector3 vtx(vertex->GetX(), vertex->GetY(), vertex->GetZ());
+
+ fH2D["vertex_y_vs_x"]->Fill(vtx[0], vtx[1]);
+ fH2D["vertex_x_vs_z"]->Fill(vtx[2], vtx[0]);
+ fH2D["vertex_y_vs_z"]->Fill(vtx[2], vtx[1]);
+
+ for (const CbmGlobalTrack* trk : fGlbTrks) {
+ float trk_tx = trk->GetParamFirst()->GetTx();
+ float trk_ty = trk->GetParamFirst()->GetTy();
+ float trk_x0 = trk->GetParamFirst()->GetX();
+ float trk_y0 = trk->GetParamFirst()->GetY();
+ float trk_z0 = trk->GetParamFirst()->GetZ();
+
+ TVector3 p(trk_x0, trk_y0, trk_z0);
+ TVector3 e(trk_tx, trk_ty, 1);
+ TVector3 trk_to_vtx = ((vtx - p).Cross(e)) * (1. / e.Mag());
+
+ fH2D["dca_dx_vs_x"]->Fill(trk_to_vtx.Px(), vtx.Px());
+ fH2D["dca_dx_vs_y"]->Fill(trk_to_vtx.Px(), vtx.Py());
+ fH2D["dca_dx_vs_z"]->Fill(trk_to_vtx.Px(), vtx.Pz());
+
+ fH2D["dca_dy_vs_x"]->Fill(trk_to_vtx.Py(), vtx.Px());
+ fH2D["dca_dy_vs_y"]->Fill(trk_to_vtx.Py(), vtx.Py());
+ fH2D["dca_dy_vs_z"]->Fill(trk_to_vtx.Py(), vtx.Pz());
+
+ fH2D["dca_dz_vs_x"]->Fill(trk_to_vtx.Pz(), vtx.Px());
+ fH2D["dca_dz_vs_y"]->Fill(trk_to_vtx.Pz(), vtx.Py());
+ fH2D["dca_dz_vs_z"]->Fill(trk_to_vtx.Pz(), vtx.Pz());
+ }
+
+ if (input_has_mc) {
+ TVector3 mc_vertex;
+ int n_of_mc_trks = fMCTrkArray->GetEntriesFast();
+ for (int mc_trk_idx = 0; mc_trk_idx < n_of_mc_trks; mc_trk_idx++) {
+ auto mc_trk = (CbmMCTrack*) fMCTrkArray->UncheckedAt(mc_trk_idx);
+ if (mc_trk->GetMotherId() == -1) {
+ mc_trk->GetStartVertex(mc_vertex);
+ break;
+ }
+ }
+
+ fH2D["res_x_vs_x"]->Fill(vtx[0] - mc_vertex[0], vtx[0]);
+ fH2D["res_x_vs_y"]->Fill(vtx[0] - mc_vertex[0], vtx[1]);
+ fH2D["res_x_vs_z"]->Fill(vtx[0] - mc_vertex[0], vtx[2]);
+
+ fH2D["res_y_vs_x"]->Fill(vtx[1] - mc_vertex[1], vtx[0]);
+ fH2D["res_y_vs_y"]->Fill(vtx[1] - mc_vertex[1], vtx[1]);
+ fH2D["res_y_vs_z"]->Fill(vtx[1] - mc_vertex[1], vtx[2]);
+
+ fH2D["res_z_vs_x"]->Fill(vtx[2] - mc_vertex[2], vtx[0]);
+ fH2D["res_z_vs_y"]->Fill(vtx[2] - mc_vertex[2], vtx[1]);
+ fH2D["res_z_vs_z"]->Fill(vtx[2] - mc_vertex[2], vtx[2]);
+ }
+ }
+
+ // Clear track containers
+ fGlbTrks.clear();
+}
+
+void CbmEventVertexDca::ProcessEvent(CbmEvent* evt)
+{
+ int nb_of_glob_trk = evt->GetNofData(ECbmDataType::kGlobalTrack);
+ fH1D["ca_track_multiplicity_all"]->Fill(nb_of_glob_trk);
+
+ if (fAnalysisCuts != nullptr && !fAnalysisCuts->CheckEvent(evt)) return;
+
+
+ for (int evt_glob_trk_idx = 0; evt_glob_trk_idx < nb_of_glob_trk; evt_glob_trk_idx++) {
+ int glob_trk_idx = evt->GetIndex(ECbmDataType::kGlobalTrack, evt_glob_trk_idx);
+ ProcessGlobalTrack((CbmGlobalTrack*) fGlbTrkArray->UncheckedAt(glob_trk_idx));
+ }
+
+ CheckVertex();
+}
+
+void CbmEventVertexDca::ProcessGlobalTrack(CbmGlobalTrack* trk)
+{
+ auto sts_trk_idx = trk->GetStsTrackIndex();
+ auto rich_trk_idx = trk->GetRichRingIndex();
+ auto much_trk_idx = trk->GetMuchTrackIndex();
+ auto trd_trk_idx = trk->GetTrdTrackIndex();
+ auto tof_trk_idx = trk->GetTofTrackIndex();
+ float trk_p_val = TMath::Prob(trk->GetChi2(), trk->GetNDF());
+
+ // Apply GlobalTracks cuts
+ int32_t mvd_trk_size = sts_trk_idx != -1 && fStsTrkArray != nullptr
+ ? ((CbmStsTrack*) fStsTrkArray->UncheckedAt(sts_trk_idx))->GetNofMvdHits()
+ : 0;
+
+ int32_t sts_trk_size = sts_trk_idx != -1 && fStsTrkArray != nullptr
+ ? ((CbmStsTrack*) fStsTrkArray->UncheckedAt(sts_trk_idx))->GetNofStsHits()
+ : 0;
+
+ int32_t rich_trk_size = rich_trk_idx != -1 && fRchTrkArray != nullptr
+ ? ((CbmTrack*) fRchTrkArray->UncheckedAt(rich_trk_idx))->GetNofHits()
+ : 0;
+
+ int32_t much_trk_size = much_trk_idx != -1 && fMchTrkArray != nullptr
+ ? ((CbmTrack*) fMchTrkArray->UncheckedAt(much_trk_idx))->GetNofHits()
+ : 0;
+
+ int32_t trd_trk_size = trd_trk_idx != -1 && fTrdTrkArray != nullptr
+ ? ((CbmTrack*) fTrdTrkArray->UncheckedAt(trd_trk_idx))->GetNofHits()
+ : 0;
+
+ int32_t tof_trk_size = tof_trk_idx != -1 && fTofTrkArray != nullptr
+ ? ((CbmTrack*) fTofTrkArray->UncheckedAt(tof_trk_idx))->GetNofHits()
+ : 0;
+
+ if (fAnalysisCuts != nullptr
+ && (!fAnalysisCuts->Check(CbmCutId::kGlobalTrackMvdSize, mvd_trk_size)
+ || !fAnalysisCuts->Check(CbmCutId::kGlobalTrackStsSize, sts_trk_size)
+ || !fAnalysisCuts->Check(CbmCutId::kGlobalTrackRichSize, rich_trk_size)
+ || !fAnalysisCuts->Check(CbmCutId::kGlobalTrackMuchSize, much_trk_size)
+ || !fAnalysisCuts->Check(CbmCutId::kGlobalTrackTrdSize, trd_trk_size)
+ || !fAnalysisCuts->Check(CbmCutId::kGlobalTrackTofSize, tof_trk_size)
+ || !fAnalysisCuts->Check(CbmCutId::kGlobalTrackChi2, trk->GetChi2())
+ || !fAnalysisCuts->Check(CbmCutId::kGlobalTrackPval, trk_p_val))) {
+ return;
+ }
+
+ fGlbTrks.push_back(trk);
+}
+
+void CbmEventVertexDca::Exec(Option_t*)
+{
+ LOG(info) << "Running CbmEventVertexDca;\t"
+ << "Entry: " << entry_;
+
+ // Check if CbmEvent
+ if (fCbmEvtArray != nullptr) {
+ uint32_t nb_events = fCbmEvtArray->GetEntriesFast();
+ LOG(info) << Form(" ... number of CbmEvent: %u\n", nb_events);
+ for (uint32_t evt_idx = 0; evt_idx < nb_events; evt_idx++) {
+ ProcessEvent((CbmEvent*) fCbmEvtArray->UncheckedAt(evt_idx));
+ }
+ }
+ else {
+
+ LOG(debug) << " ... running time-like*** Using all tracks in TS - risk of high combinatorial";
+
+ uint32_t nb_of_glob_trk = fGlbTrkArray->GetEntriesFast();
+ LOG(info) << Form("Number of GlobalTracks: %u\n", nb_of_glob_trk);
+ for (uint32_t glob_trk_idx = 0; glob_trk_idx < nb_of_glob_trk; glob_trk_idx++) {
+ CbmGlobalTrack* glob_trk = (CbmGlobalTrack*) fGlbTrkArray->UncheckedAt(glob_trk_idx);
+ ProcessGlobalTrack(glob_trk);
+ }
+
+ CheckVertex();
+ }
+ entry_++;
+}
+
+InitStatus CbmEventVertexDca::Init()
+{
+ LOG(debug) << "Init CbmEventVertexDca ...";
+
+ FairRootManager* ioman = FairRootManager::Instance();
+ if (ioman != nullptr) {
+ fCbmEvtArray = (TClonesArray*) ioman->GetObject("CbmEvent");
+
+ fGlbTrkArray = (TClonesArray*) ioman->GetObject("GlobalTrack");
+
+ fStsTrkArray = (TClonesArray*) ioman->GetObject("StsTrack");
+ fRchTrkArray = (TClonesArray*) ioman->GetObject("RichTrack");
+ fMchTrkArray = (TClonesArray*) ioman->GetObject("MuchTrack");
+ fTrdTrkArray = (TClonesArray*) ioman->GetObject("TrdTrack");
+ fTofTrkArray = (TClonesArray*) ioman->GetObject("TofTrack");
+
+ fStsHitArray = (TClonesArray*) ioman->GetObject("StsHit");
+
+ fStsCluArray = (TClonesArray*) ioman->GetObject("StsCluster");
+
+ fMCTrkArray = (TClonesArray*) ioman->GetObject("MCTrack");
+ input_has_mc = fMCTrkArray != nullptr;
+ }
+
+ LoadSetup();
+
+ BookHistograms();
+
+ fGlbTrks.reserve(100);
+
+ return kSUCCESS;
+}
+
+void CbmEventVertexDca::SetVertexFinder(std::shared_ptr<CbmDcaVertexFinder> vtx_finder) { fVertexFinder = vtx_finder; }
diff --git a/analysis/detectors/sts/CbmEventVertexDca.h b/analysis/detectors/sts/CbmEventVertexDca.h
new file mode 100644
index 0000000000000000000000000000000000000000..7ff35da8e935f248268b3374d6f0e99fc2141265
--- /dev/null
+++ b/analysis/detectors/sts/CbmEventVertexDca.h
@@ -0,0 +1,80 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+#ifndef CBMEVENTVERTEXDCA_H
+#define CBMEVENTVERTEXDCA_H
+
+#include "CbmDcaVertexFinder.h"
+#include "CbmEvent.h"
+#include "CbmGlobalTrack.h"
+#include "CbmMCTrack.h"
+#include "CbmStsAnaBase.h"
+#include "CbmStsCluster.h"
+#include "CbmStsHit.h"
+#include "CbmStsTrack.h"
+#include "CbmStsUtils.h"
+#include "CbmTrack.h"
+
+#include <FairTask.h>
+
+#include <TClonesArray.h>
+#include <TVector3.h>
+
+#include <cstring>
+#include <map>
+#include <unordered_map>
+#include <unordered_set>
+
+/** \class CbmEventVertexDca
+ * \brief Task for Event based vertex reconstruction.
+ *
+ * The current approach is to average the PCA of multiple pair of tracks.
+ * Therfore the vertex found, is a primary vertex approximation.
+ */
+class CbmEventVertexDca : public FairTask, public CbmStsAnaBase {
+ public:
+ CbmEventVertexDca() = default;
+ ~CbmEventVertexDca() = default;
+
+ InitStatus Init();
+ void Exec(Option_t*);
+ void FinishTask();
+
+ void SetVertexFinder(std::shared_ptr<CbmDcaVertexFinder>);
+
+ private:
+ std::shared_ptr<CbmDcaVertexFinder> fVertexFinder;
+
+ bool input_has_mc{false};
+
+ std::vector<CbmGlobalTrack*> fGlbTrks;
+
+ TClonesArray* fMCTrkArray{nullptr};
+ TClonesArray* fCbmEvtArray{nullptr};
+ TClonesArray* fGlbTrkArray{nullptr};
+ TClonesArray* fStsTrkArray{nullptr};
+ TClonesArray* fRchTrkArray{nullptr};
+ TClonesArray* fMchTrkArray{nullptr};
+ TClonesArray* fTrdTrkArray{nullptr};
+ TClonesArray* fTofTrkArray{nullptr};
+ TClonesArray* fStsHitArray{nullptr};
+ TClonesArray* fStsCluArray{nullptr};
+
+ void BookHistograms();
+
+ void CheckVertex();
+
+ /** \brief Process an Cbm events
+ * It filters event based on the provided CbmCutMap
+ */
+ void ProcessEvent(CbmEvent*);
+
+ /** \brief Process an Global tracks
+ * It filters the tracks based on the provided CbmCutMap
+ */
+ void ProcessGlobalTrack(CbmGlobalTrack*);
+
+ ClassDef(CbmEventVertexDca, 1);
+};
+#endif
diff --git a/analysis/detectors/sts/CbmSpillCheck.cxx b/analysis/detectors/sts/CbmSpillCheck.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..ca6f883e51522a55c5d2bf4c5bab7147a86dec61
--- /dev/null
+++ b/analysis/detectors/sts/CbmSpillCheck.cxx
@@ -0,0 +1,89 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+#include "CbmSpillCheck.h"
+
+#include <TFitResult.h>
+
+#include <ctime>
+#include <iostream>
+#include <typeinfo>
+
+CbmSpillCheck::CbmSpillCheck() : fRefSystem{ECbmModuleId::kBmon}, fRateMinPercent{0.2}, fRateMaxPercent{0.5}
+{
+ LOG(debug) << "Creating an instance of CbmSpillCheck ...";
+}
+
+CbmSpillCheck::CbmSpillCheck(ECbmModuleId ref, double lvl_min, double lvl_max)
+ : fRefSystem{ref}
+ , fRateMinPercent{lvl_min}
+ , fRateMaxPercent{lvl_max}
+{
+ LOG(debug) << "Creating an instance of CbmSpillCheck ...";
+}
+
+void CbmSpillCheck::BookHistograms(ECbmModuleId system)
+{
+ std::string h_name = Form("%s Digi Rate", ToString(system).c_str());
+ LOG(debug) << "Booking rate for : " << h_name;
+ fG1D[h_name] = std::make_unique<TGraphErrors>();
+}
+
+InitStatus CbmSpillCheck::Init()
+{
+ FairRootManager* ioman = FairRootManager::Instance();
+ if (ioman != nullptr) {
+ fDigiManager = CbmDigiManager::Instance();
+ fDigiManager->Init();
+
+ if (!fDigiManager->IsPresent(fRefSystem)) {
+ LOG(fatal) << GetName() << ": No " << ToString(fRefSystem) << " branch in input!";
+ }
+
+ for (ECbmModuleId system = ECbmModuleId::kRef; system != ECbmModuleId::kLastModule; ++system) {
+ if (fDigiManager->IsPresent(system)) BookHistograms(system);
+ }
+
+ return kSUCCESS;
+ }
+ return kERROR;
+}
+
+void CbmSpillCheck::Exec(Option_t*)
+{
+ LOG(info) << "Running CbmSpillCheck ...";
+
+ size_t nb_ref_digis = fDigiManager->GetNofDigis(fRefSystem);
+ fRateMin = fRateMin == -1 ? nb_ref_digis : std::min(nb_ref_digis, size_t(fRateMin));
+ fRateMax = fRateMax == -1 ? nb_ref_digis : std::max(nb_ref_digis, size_t(fRateMax));
+ fRefRate.push_back(nb_ref_digis);
+
+ for (ECbmModuleId system = ECbmModuleId::kRef; system != ECbmModuleId::kLastModule; ++system) {
+ if (fDigiManager->IsPresent(system)) {
+ size_t n_of_digis = fDigiManager->GetNofDigis(system);
+ std::string h_name = Form("%s Digi Rate", ToString(system).c_str());
+ fG1D[h_name]->SetPoint(fNbPoints, fNbPoints, n_of_digis);
+ }
+ }
+ fNbPoints++;
+}
+
+void CbmSpillCheck::Finish()
+{
+ double spill_lvl_off = fRateMin + fRateMinPercent * (fRateMax - fRateMin);
+ double spill_lvl_on = fRateMin + fRateMaxPercent * (fRateMax - fRateMin);
+ fSpillStatus = std::vector<int>(fRefRate.size(), 0);
+ for (size_t ts_idx = 0; ts_idx < fRefRate.size(); ts_idx++) {
+ fSpillStatus[ts_idx] = fRefRate[ts_idx] <= spill_lvl_off ? -1 : (fRefRate[ts_idx] <= spill_lvl_on ? 0 : 1);
+ std::cout << fSpillStatus[ts_idx] << std::endl;
+ }
+
+ std::cout << spill_lvl_off << std::endl;
+ std::cout << spill_lvl_on << std::endl;
+
+ SaveToFile();
+ fG1D.clear();
+ fH1D.clear();
+ fH2D.clear();
+}
diff --git a/analysis/detectors/sts/CbmSpillCheck.h b/analysis/detectors/sts/CbmSpillCheck.h
new file mode 100644
index 0000000000000000000000000000000000000000..3a93596befeadff88a63b5a997dceb03fc4e9ffe
--- /dev/null
+++ b/analysis/detectors/sts/CbmSpillCheck.h
@@ -0,0 +1,61 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+#ifndef CBMSPILLCHECK_H
+#define CBMSPILLCHECK_H
+
+#include "CbmDefs.h"
+#include "CbmDigiManager.h"
+#include "CbmStsAnaBase.h"
+
+#include <FairTask.h>
+
+#include <unordered_set>
+
+/** \class CbmSpillCheck
+ * \brief Task for checking the spill status based on a reference module.
+ *
+ * This class inherits from FairTask and CbmStsAnaBase.
+ * It provides functionality for checking the spill status based on the Digis rate of a reference module.
+ */
+class CbmSpillCheck : public FairTask, public CbmStsAnaBase {
+ public:
+ /** \brief Default constructor */
+ CbmSpillCheck();
+
+ /** \brief Parameterized constructor
+ *
+ * \param ref The ECbmModuleId representing the reference module for spill checking.
+ * \param spill_off_prcnt The percentage threshold for spill-off detection (default: 0.2).
+ * \param spill_on_prcnt The percentage threshold for spill-on detection (default: 0.5).
+ */
+ CbmSpillCheck(ECbmModuleId ref = ECbmModuleId::kBmon, double spill_off_prcnt = 0.2, double spill_on_prcnt = 0.5);
+
+ ~CbmSpillCheck() = default;
+
+ InitStatus Init();
+
+ void Exec(Option_t*);
+
+ void Finish();
+
+ private:
+ ECbmModuleId fRefSystem;
+ double fRateMinPercent;
+ double fRateMaxPercent;
+ int fRateMax{-1};
+ int fRateMin{-1};
+ int fNbPoints{0};
+ std::vector<int> fRefRate;
+ std::vector<int> fSpillStatus;
+
+
+ CbmDigiManager* fDigiManager{nullptr};
+
+ /** \brief Book histograms for a specific module */
+ void BookHistograms(ECbmModuleId);
+
+ ClassDef(CbmSpillCheck, 1);
+};
+#endif
diff --git a/analysis/detectors/sts/CbmStsAnaBase.cxx b/analysis/detectors/sts/CbmStsAnaBase.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..3076f5af47bdf5953681d378027e0fe7ed3d2253
--- /dev/null
+++ b/analysis/detectors/sts/CbmStsAnaBase.cxx
@@ -0,0 +1,123 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+#include "CbmStsAnaBase.h"
+
+#include "CbmStsModule.h"
+#include "CbmStsStation.h"
+#include "FairRootManager.h"
+
+#include <TGeoBBox.h> // for Geometry loading from file
+#include <TGeoManager.h> // for Geometry loading from file
+#include <TGeoMatrix.h> // for Geometry loading from file
+#include <TGeoNode.h> // for Geometry loading from file
+#include <TGeoPhysicalNode.h> // for Geometry loading from file
+#include <TGeoVolume.h> // for Geometry loading from file
+
+#include <cassert>
+
+
+void CbmStsAnaBase::SetCutMap(CbmCutMap* cuts_map)
+{
+ assert(cuts_map != nullptr);
+ fAnalysisCuts = cuts_map;
+}
+
+void CbmStsAnaBase::LoadSetup()
+{
+ CbmStsSetup* sts_setup = CbmStsSetup::Instance();
+ if (!sts_setup->IsInit()) {
+ sts_setup->Init(nullptr);
+ }
+
+ if (sts_setup == nullptr) {
+ LOG(warning) << "Something is wrong with the setup ...";
+ return;
+ }
+
+ nb_sts_station_ = sts_setup->GetNofStations();
+ assert(nb_sts_station_ != 0);
+
+ for (unsigned short unit_idx = 0; unit_idx < nb_sts_station_; unit_idx++) {
+ CbmStsStation* sts_unit = sts_setup->GetStation(unit_idx);
+ fStsGeoInfo[unit_idx] = std::vector<double>(5, 0);
+ fStsGeoInfo[unit_idx][0] = sts_unit->GetXmin();
+ fStsGeoInfo[unit_idx][1] = sts_unit->GetXmax();
+ fStsGeoInfo[unit_idx][2] = sts_unit->GetYmin();
+ fStsGeoInfo[unit_idx][3] = sts_unit->GetYmax();
+ fStsGeoInfo[unit_idx][4] = sts_unit->GetZ();
+ }
+
+ for (unsigned short module_idx = 0; module_idx < sts_setup->GetNofModules(); module_idx++) { // Loop over modules
+ CbmStsModule* sts_module = sts_setup->GetModule(module_idx);
+ assert(sts_module);
+ int32_t address = sts_module->GetAddress();
+ assert(sts_module->GetNofDaughters() == 1);
+
+ CbmStsElement* sts_sensor = sts_module->GetDaughter(0);
+ TGeoPhysicalNode* sensor_node = sts_sensor->GetPnode();
+ TGeoVolume* sensor_volume = sensor_node->GetVolume();
+ TGeoBBox* sensor_shape = (TGeoBBox*) sensor_volume->GetShape();
+ TGeoMatrix* sensor_matrix = sensor_node->GetMatrix();
+
+ // const double* local = sensor_shape->GetOrigin();
+ const double* trans = sensor_matrix->GetTranslation();
+ fStsGeoInfo[address] = std::vector<double>(6, 0);
+ fStsGeoInfo[address][0] = trans[0] - sensor_shape->GetDX(); // Xmin
+ fStsGeoInfo[address][1] = trans[0] + sensor_shape->GetDX(); // Xmax
+ fStsGeoInfo[address][2] = trans[1] - sensor_shape->GetDY(); // Ymin
+ fStsGeoInfo[address][3] = trans[1] + sensor_shape->GetDY(); // Ymax
+ fStsGeoInfo[address][4] = trans[2] - sensor_shape->GetDZ(); // Zmin
+ fStsGeoInfo[address][5] = trans[2] + sensor_shape->GetDZ(); // Zmax
+
+ if (address > 8) {
+ switch (2 * (int) sensor_shape->GetDY()) {
+ case 2: fFirstZStrip[address] = 2048 - 32; break;
+ case 4: fFirstZStrip[address] = 2048 - 88; break;
+ case 6: fFirstZStrip[address] = 2048 - 134; break;
+ case 12: fFirstZStrip[address] = 2048 - 274; break;
+ default:
+ LOG(warning) << Form("Unknown sensor shape 0x%x: [%0.3f , %0.3f , %0.3f]", address, sensor_shape->GetDX(),
+ sensor_shape->GetDY(), sensor_shape->GetDZ());
+ break;
+ }
+ }
+ }
+
+ LOG(debug) << "Loaded setup information:\n";
+ for (auto& [address, geo] : fStsGeoInfo) {
+ LOG(debug) << Form("0x%x: [%+0.3f, %+0.3f] ^ [%+0.3f, %+0.3f] ^ [%+0.3f, %+0.3f]", address, geo[0], geo[1], geo[2],
+ geo[3], geo[4], geo[5]);
+ }
+}
+
+void CbmStsAnaBase::SaveToFile()
+{
+ FairRootManager* ioman = FairRootManager::Instance();
+ auto sink = ioman->GetSink();
+
+ if (!sink) {
+ LOG(fatal) << "Check the configuration: sink is invalid!";
+ return;
+ }
+
+ LOG(info) << Form("Task histograms will be saved to %s ...", sink->GetFileName().Data());
+
+ for (auto& [name, gr] : fG1D) {
+ if (gr == nullptr) continue;
+ sink->WriteObject(gr.get(), name.c_str());
+ }
+ for (auto& [name, h] : fH1D) {
+ if (h == nullptr) continue;
+ sink->WriteObject(h.get(), name.c_str());
+ }
+ for (auto& [name, h] : fH2D) {
+ if (h == nullptr) continue;
+ sink->WriteObject(h.get(), name.c_str());
+ }
+ for (auto& [name, h] : fH2DShared) {
+ if (h == nullptr) continue;
+ sink->WriteObject(h.get(), name.c_str());
+ }
+}
diff --git a/analysis/detectors/sts/CbmStsAnaBase.h b/analysis/detectors/sts/CbmStsAnaBase.h
new file mode 100644
index 0000000000000000000000000000000000000000..1dfd3941fc97780efc98d66bd6bdd205da7f3c95
--- /dev/null
+++ b/analysis/detectors/sts/CbmStsAnaBase.h
@@ -0,0 +1,101 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+#ifndef CBMSTSANABASE_H
+#define CBMSTSANABASE_H
+
+#include "CbmCutMap.h"
+#include "CbmDefs.h"
+#include "CbmStsSetup.h"
+#include "CbmStsUtils.h"
+
+#include <FairRootManager.h>
+#include <Logger.h>
+
+#include <TCanvas.h>
+#include <TClonesArray.h>
+#include <TFile.h>
+#include <TGraphErrors.h>
+#include <TH1D.h>
+#include <TH2D.h>
+#include <TSystem.h>
+
+#include <cstring>
+#include <map>
+#include <optional>
+#include <unordered_set>
+
+
+enum InputType
+{
+ kMCEventMode,
+ kMCTimeMode,
+ kEventMode,
+ kTimeMode
+};
+
+/**
+ * @class CbmStsAnaBase
+ * @brief Base class for STS analysis.
+ *
+ * This class provides a base for STS (Silicon Tracking System) analysis. It contains
+ * common functionalities and data members used by derived analysis classes.
+ */
+class CbmStsAnaBase {
+ public:
+ CbmStsAnaBase() = default;
+ ~CbmStsAnaBase() = default;
+
+ /**
+ * @brief Set the cut map for analysis.
+ * @param cutMap Pointer to the CbmCutMap object.
+ */
+ void SetCutMap(CbmCutMap* map);
+
+ /**
+ * @brief User defined sensor translations.
+ * @param user_mat Input translations.
+ */
+ void UserAlignment(const std::map<int32_t, std::vector<double>>& user_mat) { fUserAlignment = user_mat; }
+
+ /**
+ * @brief Virtual function to draw analysis results.
+ */
+ virtual void DrawResults() {}
+
+ protected:
+ uint entry_{0};
+ std::unique_ptr<TFile> fReportFile;
+ std::unordered_set<int32_t> fAddressBook;
+ std::map<std::string, std::unique_ptr<TGraphErrors>> fG1D; // Map of TGraphErrors objects.
+ std::map<std::string, std::unique_ptr<TH1D>> fH1D; // Map of TH1D objects.
+ std::map<std::string, std::unique_ptr<TH2D>> fH2D; // Map of TH2D objects.
+ std::map<std::string, std::shared_ptr<TH2D>> fH2DShared; // Map of TH2D objects.
+ std::map<std::string, std::unique_ptr<TCanvas>> fCanvas; // Map of TH2D objects.
+
+ int nb_sts_station_{8}; // Number of STS stations.
+ std::unordered_map<int32_t, std::vector<double>> fStsGeoInfo; // Map of STS geometry information.
+
+ std::map<int32_t, std::vector<double>> fUserAlignment; // Stores the user-defined alignment information.
+
+ std::unordered_map<int32_t, int> fFirstZStrip;
+
+ CbmCutMap* fAnalysisCuts{nullptr};
+
+ int fRunId{-1};
+
+ /**
+ * @brief Load the STS setup and fill the map with XYZ boundaries for each STS setup element.
+ * It maps the first z strip for each module depending on the size of the sensor.
+ */
+ void LoadSetup();
+
+ /**
+ * @brief It write all mapped objects to the FairRunAna sink file
+ */
+ void SaveToFile();
+
+ ClassDef(CbmStsAnaBase, 1);
+};
+#endif
diff --git a/analysis/detectors/sts/CbmStsAnaLinkDef.h b/analysis/detectors/sts/CbmStsAnaLinkDef.h
index 6b334db83021bf24d9e82831896d7088df2a815e..8f7a841e294a4c45eca0e17b9eba9634d1f06471 100644
--- a/analysis/detectors/sts/CbmStsAnaLinkDef.h
+++ b/analysis/detectors/sts/CbmStsAnaLinkDef.h
@@ -10,5 +10,20 @@
#pragma link C++ class CbmStsWkn + ;
+#pragma link C++ class CbmCutMap + ;
+#pragma link C++ class CbmStsAnaBase + ;
+#pragma link C++ namespace CbmStsUtils + ;
+#pragma link C++ class CbmStsAnalysis + ;
+#pragma link C++ class CbmSpillCheck + ;
+#pragma link C++ class CbmStsChannelQA + ;
+#pragma link C++ class CbmStsTimeCal + ;
+#pragma link C++ class CbmStsHitAna + ;
+#pragma link C++ class CbmStsCorrelation + ;
+#pragma link C++ class CbmStsRecoBeamSpot + ;
+#pragma link C++ class CbmDcaVertexFinder + ;
+#pragma link C++ class CbmStsEfficiency + ;
+#pragma link C++ class CbmStsResolution + ;
+#pragma link C++ class CbmEventVertexDca + ;
+
#endif /* __CINT__ */
diff --git a/analysis/detectors/sts/CbmStsAnalysis.cxx b/analysis/detectors/sts/CbmStsAnalysis.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..497afb8ee6ba0bf3160efd17560b47dbaf3379d1
--- /dev/null
+++ b/analysis/detectors/sts/CbmStsAnalysis.cxx
@@ -0,0 +1,12 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+#include "CbmStsAnalysis.h"
+
+#include "TList.h"
+
+#include <cassert>
+#include <iomanip>
+
+void CbmStsAnalysis::SetAnalysisCuts(CbmCutMap* cuts_map) { fAnaCutsMap = cuts_map; }
diff --git a/analysis/detectors/sts/CbmStsAnalysis.h b/analysis/detectors/sts/CbmStsAnalysis.h
new file mode 100644
index 0000000000000000000000000000000000000000..f5d7b2e259fef1bd676025d2307235209f885d8b
--- /dev/null
+++ b/analysis/detectors/sts/CbmStsAnalysis.h
@@ -0,0 +1,39 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+#ifndef CBMSTSANALYSIS_H
+#define CBMSTSANALYSIS_H
+
+// Includes from C/C++
+#include <unordered_map>
+#include <unordered_set>
+
+// Includes from FairRoot
+#include <FairRunAna.h>
+#include <Logger.h>
+
+// Includes from CbmRoot
+#include "CbmStsAnaBase.h"
+#include "CbmStsDigi.h"
+
+// Includes from ROOT
+#include <TClonesArray.h>
+#include <TFile.h>
+#include <TSystem.h>
+#include <TTree.h>
+
+class CbmStsAnalysis : public FairRunAna {
+ public:
+ CbmStsAnalysis() = default;
+ ~CbmStsAnalysis() = default;
+
+ void SetAnalysisCuts(CbmCutMap*);
+
+ private:
+ CbmCutMap* fAnaCutsMap{nullptr};
+
+ ClassDef(CbmStsAnalysis, 1);
+};
+
+#endif
diff --git a/analysis/detectors/sts/CbmStsChannelQA.cxx b/analysis/detectors/sts/CbmStsChannelQA.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..d2734a216db6210090e9ec8e713d0f8a20601f89
--- /dev/null
+++ b/analysis/detectors/sts/CbmStsChannelQA.cxx
@@ -0,0 +1,315 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+#include "CbmStsChannelQA.h"
+
+#include "CbmStsDigi.h"
+#include "TArrow.h"
+#include "TBox.h"
+#include "TPaveLabel.h"
+#include "TPaveText.h"
+#include "TText.h"
+
+#include <TParameter.h>
+
+#include <ctime>
+#include <iostream>
+#include <typeinfo>
+
+CbmStsChannelQA::CbmStsChannelQA(int dead_threshold) : fActiveMinEntries(dead_threshold) {}
+CbmStsChannelQA::CbmStsChannelQA(int dead_threshold, double noise_threshold,
+ std::optional<std::pair<size_t, size_t>> spill_on_off_threshold)
+ : fActiveMinEntries(dead_threshold)
+ , fSBThreshold(noise_threshold)
+ , fSpillThresholds(spill_on_off_threshold)
+{
+ assert(noise_threshold > 0);
+ if (fSpillThresholds.has_value()) {
+ assert(fSpillThresholds->first <= fSpillThresholds->second);
+ fSpillSections = {":ramp", ":spill_on", ":spill_off"};
+ }
+}
+
+InitStatus CbmStsChannelQA::Init()
+{
+ FairRootManager* ioman = FairRootManager::Instance();
+ if (ioman != nullptr) {
+ fDigiManager = CbmDigiManager::Instance();
+ fDigiManager->Init();
+
+ if (!fDigiManager->IsPresent(ECbmModuleId::kSts)) LOG(fatal) << GetName() << ": No StsDigi branch in input!";
+
+ LOG(info) << "CbmStsChannelQA configuration:";
+ LOG(info) << "Report lvl: " << fReportLvl;
+ LOG(info) << "Active channel minimum entries: " << fActiveMinEntries;
+ LOG(info) << "Noisy channel S/B thresholds: " << fSBThreshold;
+ if (fSpillThresholds.has_value()) {
+ LOG(info) << Form("Spill OFF: RefDet_DigiRate < %lu", fSpillThresholds->first);
+ LOG(info) << Form("Spill ON: RefDet_DigiRate > %lu", fSpillThresholds->second);
+ }
+ else {
+ LOG(info) << "No Spill section defined. Using all time slices";
+ }
+
+ return kSUCCESS;
+ }
+ return kERROR;
+}
+
+
+void CbmStsChannelQA::BookHistograms(int32_t address)
+{
+ LOG(debug) << Form("Booking histograms for module_addr: 0x%x", address);
+
+ double dt_max = 120; // HARDCODE
+ auto t_binning = cbm_sts_utils::HBinning{uint32_t(dt_max / 3.125), 0, dt_max};
+
+ std::string h_name;
+ for (auto modifier : fSpillSections) {
+ h_name = Form("0x%x_charge_vs_channel%s", address, modifier);
+ LOG(debug) << h_name;
+ fH2D[h_name] = std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), 2048, 0, 2048, 31, 1, 32);
+ fH2D[h_name]->GetXaxis()->SetTitle("Channel_{StsDigi}");
+ fH2D[h_name]->GetYaxis()->SetTitle("Charge_{StsDigi} [ADC]");
+
+ LOG(debug) << h_name;
+ h_name = Form("0x%x_dt_vs_charge%s", address, modifier);
+ fH2D[h_name] = std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), 31, 1, 32, t_binning.n_of_bins,
+ t_binning.x_min, t_binning.x_max);
+ fH2D[h_name]->GetYaxis()->SetTitle("Charge_{StsDigi} [ADC]");
+ fH2D[h_name]->GetYaxis()->SetTitle("Time difference [ns]");
+
+ LOG(debug) << h_name;
+ h_name = Form("0x%x_dt_vs_channel%s", address, modifier);
+ fH2D[h_name] = std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), 2048, 0, 2048, t_binning.n_of_bins,
+ t_binning.x_min, t_binning.x_max);
+ fH2D[h_name]->GetXaxis()->SetTitle("Channel_{StsDigi}");
+ fH2D[h_name]->GetYaxis()->SetTitle("Time difference [ns]");
+ }
+
+ fAddressBook.insert(address);
+}
+
+void CbmStsChannelQA::ResetLastTime()
+{
+ for (auto& [address, time] : fLastDigiTime) {
+ for (int chn = 0; chn < 2048; chn++) {
+ time[chn] = -1;
+ }
+ }
+}
+
+void CbmStsChannelQA::Exec(Option_t*)
+{
+ LOG(info) << "Running CbmStsChannelQA ...";
+
+ // Get Spill status
+ size_t nb_ref_digis = fDigiManager->GetNofDigis(ECbmModuleId::kBmon);
+
+ /**
+ * -1 Off spill
+ * 0 Spill ramp
+ * +1 On spill
+ **/
+ int fSpillStatus = 0;
+ if (fSpillThresholds.has_value()) {
+ fSpillStatus = nb_ref_digis <= fSpillThresholds->first ? -1 : (nb_ref_digis <= fSpillThresholds->second ? 0 : 1);
+ }
+
+ std::string str_spill;
+ switch (fSpillStatus) {
+ case +1: {
+ fNbTsSpillOn++;
+ str_spill = fSpillSections[1];
+ break;
+ }
+ case -1: {
+ fNbTsSpillOff++;
+ str_spill = fSpillSections[2];
+ break;
+ }
+ default: {
+ str_spill = fSpillSections[0];
+ break;
+ }
+ }
+
+ LOG(info) << Form("TS %d\t-\t Spill section: %s", entry_, str_spill.c_str());
+
+
+ auto sts_digis_ = fDigiManager->GetArray<CbmStsDigi>();
+ size_t nb_sts_digis = fDigiManager->GetNofDigis(ECbmModuleId::kSts);
+ for (size_t sts_digi_idx = 0; sts_digi_idx < nb_sts_digis; sts_digi_idx++) { // sts digi loop
+ const CbmStsDigi* sts_digi = &sts_digis_[sts_digi_idx];
+ int32_t sts_digi_addr = sts_digi->GetAddress();
+ int32_t sts_digi_chan = sts_digi->GetChannel();
+ int32_t sts_digi_char = sts_digi->GetCharge();
+ double sts_digi_time = sts_digi->GetTime();
+
+ if (!fAddressBook.count(sts_digi_addr)) {
+ BookHistograms(sts_digi_addr);
+
+ for (int chn = 0; chn < 2048; chn++) {
+ fLastDigiTime[sts_digi_addr][chn] = -1;
+ }
+ }
+
+ fH2D[Form("0x%x_charge_vs_channel%s", sts_digi_addr, str_spill.c_str())]->Fill(sts_digi_chan, sts_digi_char);
+
+ if (fLastDigiTime[sts_digi_addr][sts_digi_chan] > 0) {
+ fH2D[Form("0x%x_dt_vs_charge%s", sts_digi_addr, str_spill.c_str())]->Fill(
+ sts_digi_char, sts_digi_time - fLastDigiTime[sts_digi_addr][sts_digi_chan]);
+ fH2D[Form("0x%x_dt_vs_channel%s", sts_digi_addr, str_spill.c_str())]->Fill(
+ sts_digi_chan, sts_digi_time - fLastDigiTime[sts_digi_addr][sts_digi_chan]);
+ }
+
+ fLastDigiTime[sts_digi_addr][sts_digi_chan] = sts_digi_time;
+
+ } // end sts digi loop
+
+ ResetLastTime();
+ entry_++;
+}
+
+void CbmStsChannelQA::CheckDeadChannels()
+{
+ const char* mod = fSpillThresholds.has_value() ? fSpillSections[1] : fSpillSections[0];
+ for (auto& module_addr : fAddressBook) {
+ auto h = (TH1D*) fH2D[Form("0x%x_charge_vs_channel%s", module_addr, mod)]->ProjectionX();
+ for (int chn = 0; chn < 2048; chn++) {
+ double entries = h->GetBinContent(chn + 1);
+
+ if (entries < fActiveMinEntries) {
+ fDeadChannelList[module_addr].push_back(chn);
+ }
+ }
+ }
+}
+
+void CbmStsChannelQA::CheckNoisyChannels()
+{
+ for (auto& module_addr : fAddressBook) {
+ TH1D* h_sgn = fH2D[Form("0x%x_charge_vs_channel:spill_on", module_addr)]->ProjectionX();
+ TH1D* h_bkg = fH2D[Form("0x%x_charge_vs_channel:spill_off", module_addr)]->ProjectionX();
+
+ h_sgn->Scale(1. / fNbTsSpillOn);
+ h_bkg->Scale(1. / fNbTsSpillOff);
+ h_sgn->Add(h_bkg, -1);
+
+ for (int chn = 0; chn < 2048; chn++) {
+ double sgn = h_sgn->GetBinContent(chn + 1);
+ double bkg = h_bkg->GetBinContent(chn + 1);
+
+ if (bkg == 0) continue;
+
+ double sb_ratio = sgn / bkg;
+
+ if (sb_ratio < fSBThreshold) {
+ LOG(debug) << Form("Noisy channel 0x%x: %d, \tS/B: %0.4f", module_addr, chn, sb_ratio);
+ fNoisyChannelList[module_addr].push_back(chn);
+ }
+ }
+ }
+}
+
+void CbmStsChannelQA::GenerateReport()
+{
+ std::string f_name = "cbm_sts_channel_qa_report.root";
+ std::unique_ptr<TFile> o_root_file = std::make_unique<TFile>(f_name.c_str(), "RECREATE");
+ int pad_size_x = 1000;
+ int pad_size_y = 1000;
+ for (auto& module_addr : fAddressBook) {
+ std::string c_name = Form("0x%x_channel", module_addr);
+ std::unique_ptr<TCanvas> c = std::make_unique<TCanvas>(c_name.c_str(), c_name.c_str(), pad_size_x, pad_size_y);
+ TH1D* h_sgn = fH2D[Form("0x%x_charge_vs_channel:spill_on", module_addr)]->ProjectionX();
+ TH1D* h_bkg = fH2D[Form("0x%x_charge_vs_channel:spill_off", module_addr)]->ProjectionX();
+
+ LOG(debug) << Form("Building report for 0x%x", module_addr);
+ h_sgn->Scale(1. / fNbTsSpillOn);
+ h_bkg->Scale(1. / fNbTsSpillOff);
+
+ h_sgn->SetLineColor(kRed);
+ h_bkg->SetLineColor(kBlack);
+
+ h_sgn->SetTitle("");
+ h_bkg->SetTitle("");
+
+ h_sgn->SetFillColorAlpha(kRed, 0.2);
+ h_bkg->SetFillColorAlpha(kBlack, 0.2);
+
+ double h_sgn_max = h_sgn->GetMaximum();
+ double h_bkg_max = h_bkg->GetMaximum();
+ double max = std::max(h_sgn_max, h_bkg_max);
+
+ h_sgn->SetMaximum(1.05 * max);
+ h_bkg->SetMaximum(1.05 * max);
+
+ c->cd(1);
+ gPad->SetLogy(1);
+ gPad->SetLeftMargin(0.15);
+ gPad->SetBottomMargin(0.12);
+ gPad->SetTopMargin(0.10);
+ gPad->SetRightMargin(0.12);
+
+ h_sgn->Draw("histo");
+ h_bkg->Draw("histo same");
+
+ // Draw a line for bad channels
+ if (fReportLvl > 1) {
+ LOG(debug) << Form("Drawing %d noisy channel markers: %lu", module_addr, fNoisyChannelList[module_addr].size());
+ for (auto& chn : fNoisyChannelList[module_addr]) {
+ TBox chn_box = TBox(chn + 0.2, 0.99 * max, chn + 0.8, 1.01 * max);
+ chn_box.SetLineColor(kBlue);
+ chn_box.SetFillColor(kBlue);
+ chn_box.DrawClone();
+ }
+ }
+
+ c->Write(Form("0x%x_channel:spill_on_vs_off.png", module_addr));
+ }
+}
+
+void CbmStsChannelQA::Finish()
+{
+ CheckDeadChannels();
+ if (fSpillThresholds.has_value()) CheckNoisyChannels();
+
+ if (fReportLvl) GenerateReport();
+
+ SaveToFile();
+
+ fH1D.clear();
+ fH2D.clear();
+
+ // Write normalization factor to file
+ std::unique_ptr<TFile> out_file = std::make_unique<TFile>("cbm_sts_channel_qa.root", "UPDATE");
+ TParameter<int>("fNbTsSpillOn", fNbTsSpillOn).Write();
+ TParameter<int>("fNbTsSpillOff", fNbTsSpillOff).Write();
+
+ std::ofstream dead_info("sts_dead_channels.info");
+ std::ofstream dead_par("sts_dead_channels.par");
+ dead_info << Form("# Module\t Number of dead channels\n");
+ int n_of_dead_chn = 0;
+ for (auto& [module_addr, list] : fDeadChannelList) {
+ dead_info << Form("0x%x\t%lu\n", module_addr, list.size());
+ for (int& dead_chn : list) {
+ dead_par << module_addr << "\t" << dead_chn << std::endl;
+ }
+ n_of_dead_chn += list.size();
+ }
+ dead_info << Form("# Total Number of dead channels: %d\n", n_of_dead_chn);
+
+ std::ofstream noisy_info("sts_noisy_channels.info");
+ std::ofstream noisy_par("sts_noisy_channels.par");
+ noisy_info << Form("# Module\t Number of noisy channels\n");
+ int n_of_noisy_chn = 0;
+ for (auto& [module_addr, list] : fNoisyChannelList) {
+ noisy_info << Form("0x%x\t%lu\n", module_addr, list.size());
+ for (int& noisy_chn : list) {
+ noisy_par << module_addr << "\t" << noisy_chn << std::endl;
+ }
+ n_of_noisy_chn += list.size();
+ }
+ noisy_info << Form("# Total Number of noisy channels: %d\n", n_of_noisy_chn);
+}
diff --git a/analysis/detectors/sts/CbmStsChannelQA.h b/analysis/detectors/sts/CbmStsChannelQA.h
new file mode 100644
index 0000000000000000000000000000000000000000..d947d5fb567894d4c018beac174940b0b3716c21
--- /dev/null
+++ b/analysis/detectors/sts/CbmStsChannelQA.h
@@ -0,0 +1,136 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+
+/**
+ * @class CbmStsChannelQA
+ * @brief FairTask for StsChannelQa.
+ *
+ * It processes StsDigis to produce two sets of output files:
+ * - sts_dead_channels.info: Summary of the amount of dead channels for each STS module (.par file with two columns <module_address> <dead_channel>).
+ * - sts_noisy_channels.info: Summary of the amount of noisy channels for each STS module (.par file with two columns <module_address> <noisy_channel>).
+ *
+ * Moreover, all the histograms/graphs generated will be dumped to: cbm_sts_channel_qa.root.
+ */
+
+#ifndef CBMSTSCHANNELQA_H
+#define CBMSTSCHANNELQA_H
+
+#include "CbmDigiManager.h"
+#include "CbmStsAnaBase.h"
+
+#include <FairTask.h>
+
+#include <TClonesArray.h>
+
+#include <cstring>
+#include <map>
+#include <unordered_set>
+
+/**
+ * @class CbmStsChannelQA
+ * @brief Quality assurance task for STS channel monitoring.
+ *
+ * This class inherits from FairTask and CbmStsAnaBase, and provides tools
+ * to analyze and monitor the behavior of STS channels during data processing.
+ *
+ * The main functions of this class include:
+ *
+ * - **Dead Channel Detection**:
+ * Identifies channels with an entry count below a minimum threshold,
+ * classifying them as non-functional or "dead".
+ *
+ * - **Noisy Channel Detection**:
+ * Evaluates the signal-to-background (S/B) ratio using the ratio of
+ * entries in-spill versus off-spill. Channels with an user-defined
+ * S/B ratio are marked as noisy.
+ *
+ * - **Charge vs. Channel Histograms**:
+ * Generates 2D histograms of charge versus channel, split by spill sections (see CbmSpillCheck documentation):
+ * - All events : if not spill information is provided
+ * - Ramp section: between the lower and upper thresholds
+ * - In-spill section : bellow the lower threshold
+ * - Off-spill section : above the upper threshold
+ *
+ * If spill information is provided the number of in|off spill TS analyzed are save to enable normalization
+ *
+ * Output files are produced:
+ * - **Summary files**:
+ * - sts_dead_channels.info: Summary of the amount of dead channels for each STS module (.par file with two columns <module_address> <dead_channel>).
+ * - sts_noisy_channels.info: Summary of the amount of noisy channels for each STS module (.par file with two columns <module_address> <noisy_channel>).
+ * - **Extendend information**:
+ * - cbm_sts_channel_qa.root: Contains all the histograms/graphs generated.
+ * - sts_dead_channels.par: Contains the list of dead channels for each module.
+ * - sts_noisy_channels.par: Contains the list of noisy channels for each module.
+ */
+class CbmStsChannelQA : public FairTask, public CbmStsAnaBase {
+ public:
+ CbmStsChannelQA() = default;
+
+ /**
+ * @brief Constructor with threshold for inactive channels.
+ * @param dead_threshold Minimum amount of entries to consider a channel active.
+ * Useful to generate the channel masking map for MC input
+ */
+ CbmStsChannelQA(int dead_threshold);
+
+ /**
+ * @brief Constructor with thresholds for dead and noisy channels. Additionally set thresholds
+ * to define OFF | Rise-Fall | ON spill status base on RefDet digi rate
+ * @param dead_threshold Minimum amount of entries to consider a channel alive.
+ * @param noise_threshold Minimum value for Signal/Background ratio to consider a channel to be noisy
+ * @param spill_on_off_threshold (Minimum,Maximum) value for RefDet digi rate to consider spill (ON,OFF) [optional analysis branching]
+ */
+ CbmStsChannelQA(int dead_threshold, double noise_threshold, const std::optional<std::pair<size_t, size_t>>);
+
+ ~CbmStsChannelQA() = default;
+
+ InitStatus Init();
+ void Exec(Option_t*);
+ void Finish();
+
+ void SetReportLevel(int lvl) { fReportLvl = lvl; }
+
+ private:
+ int fReportLvl{0};
+ const int fActiveMinEntries{1};
+ const double fSBThreshold{0.5};
+ const std::optional<std::pair<size_t, size_t>> fSpillThresholds{std::nullopt};
+
+ std::vector<const char*> fSpillSections = {":all"};
+
+ int fNbTsSpillOn{0};
+ int fNbTsSpillOff{0};
+
+ std::map<int32_t, std::vector<int>> fDeadChannelList;
+ std::map<int32_t, std::vector<int>> fNoisyChannelList;
+
+ std::map<int32_t, double[2048]> fLastDigiTime;
+ void ResetLastTime();
+
+ CbmDigiManager* fDigiManager{nullptr};
+
+ /**
+ * @brief Book histograms.
+ * @param address The module address for which histograms will be booked.
+ */
+ void BookHistograms(int32_t);
+
+ /**
+ * @brief Check for dead channels.
+ * Fills fDeadChannelList map
+ */
+ void CheckDeadChannels();
+
+ /**
+ * @brief Check for dead channels.
+ * Fills fNoisyChannelList map
+ */
+ void CheckNoisyChannels();
+
+ void GenerateReport();
+
+ ClassDef(CbmStsChannelQA, 1);
+};
+#endif
diff --git a/analysis/detectors/sts/CbmStsCorrelation.cxx b/analysis/detectors/sts/CbmStsCorrelation.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..9f0e34f32525d46f245753193da8fe7a10689cdd
--- /dev/null
+++ b/analysis/detectors/sts/CbmStsCorrelation.cxx
@@ -0,0 +1,205 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+#include "CbmStsCorrelation.h"
+
+#include "CbmPixelHit.h"
+#include "CbmStsAddress.h"
+#include "TPaveStats.h"
+#include "TStyle.h"
+
+void CbmStsCorrelation::BookHistograms()
+{
+
+ for (auto& [element_i, geo_i] : fStsGeoInfo) {
+ for (auto& [element_j, geo_j] : fStsGeoInfo) {
+
+ std::string h_name_base = "";
+ std::string obj_x, obj_y;
+ if (element_i < 8 && element_j < 8) {
+ if (element_i >= element_j) continue;
+ obj_x = Form("Sts%d", element_i);
+ obj_y = Form("Sts%d", element_j);
+ h_name_base = Form("Sts%d:Sts%d", element_i, element_j);
+ }
+ else if (element_i > 8 && element_j > 8) {
+ int32_t unit_i = CbmStsAddress::GetElementId(element_i, kStsUnit);
+ int32_t unit_j = CbmStsAddress::GetElementId(element_j, kStsUnit);
+
+ if (unit_i >= unit_j) continue;
+
+ obj_x = Form("Sts0x%x", element_i);
+ obj_y = Form("Sts0x%x", element_j);
+ h_name_base = Form("Sts0x%x:Sts0x%x", element_i, element_j);
+ }
+
+ if (!h_name_base.length()) continue;
+ LOG(debug) << Form("Booking for %s", h_name_base.c_str());
+
+ double x_min_i = geo_i[0];
+ double x_max_i = geo_i[1];
+ double y_min_i = geo_i[2];
+ double y_max_i = geo_i[3];
+
+ unsigned int nb_bins_x_i = (x_max_i - x_min_i) / cbm_sts_utils::kStsDx;
+ unsigned int nb_bins_y_i = (y_max_i - y_min_i) / cbm_sts_utils::kStsDy;
+
+ double x_min_j = geo_j[0];
+ double x_max_j = geo_j[1];
+ double y_min_j = geo_j[2];
+ double y_max_j = geo_j[3];
+
+ LOG(debug) << Form("XX_[%0.3f, %0.3f]:[%0.3f, %0.3f]", x_min_i, x_max_i, x_min_j, x_max_j);
+ LOG(debug) << Form("YY_[%0.3f, %0.3f]:[%0.3f, %0.3f]", y_min_i, y_max_i, y_min_j, y_max_j);
+
+ unsigned int nb_bins_x_j = (x_max_j - x_min_j) / cbm_sts_utils::kStsDx;
+ unsigned int nb_bins_y_j = (y_max_j - y_min_j) / cbm_sts_utils::kStsDy;
+
+ std::string h_name;
+ h_name = Form("%s_XX", h_name_base.c_str());
+ fH2D[h_name] =
+ std::make_unique<TH2D>(h_name.c_str(), "", nb_bins_x_i, x_min_i, x_max_i, nb_bins_x_j, x_min_j, x_max_j);
+ fH2D[h_name]->GetXaxis()->SetTitle(Form("X_{%s} [cm]", obj_x.c_str()));
+ fH2D[h_name]->GetYaxis()->SetTitle(Form("X_{%s} [cm]", obj_y.c_str()));
+ fH2D[h_name]->GetXaxis()->CenterTitle(true);
+ fH2D[h_name]->GetYaxis()->CenterTitle(true);
+
+ h_name = Form("%s_YY", h_name_base.c_str());
+ fH2D[h_name] =
+ std::make_unique<TH2D>(h_name.c_str(), "", nb_bins_y_i, y_min_i, y_max_i, nb_bins_y_j, y_min_j, y_max_j);
+ fH2D[h_name]->GetXaxis()->SetTitle(Form("Y_{%s} [cm]", obj_x.c_str()));
+ fH2D[h_name]->GetYaxis()->SetTitle(Form("Y_{%s} [cm]", obj_y.c_str()));
+ fH2D[h_name]->GetXaxis()->CenterTitle(true);
+ fH2D[h_name]->GetYaxis()->CenterTitle(true);
+ }
+ }
+}
+
+void CbmStsCorrelation::BuildCorrelation()
+{
+ bool is_mc_data = false;
+ if (!fStsHits.size()) return;
+ for (int sts_unit_i = 0; sts_unit_i < nb_sts_station_ - 1; sts_unit_i++) {
+ for (int sts_unit_j = sts_unit_i + 1; sts_unit_j < nb_sts_station_; sts_unit_j++) {
+ LOG(debug) << Form("Building correlations for STS%d:STS%d: %ld\t%ld", sts_unit_i, sts_unit_j,
+ fStsHits[sts_unit_i].size(), fStsHits[sts_unit_j].size());
+ for (auto sts_hit_i : fStsHits[sts_unit_i]) {
+ int32_t sts_address_i = sts_hit_i->GetAddress();
+ double sts_i_x = sts_hit_i->GetX();
+ double sts_i_y = sts_hit_i->GetY();
+
+ CbmStsHit* closest_hit = nullptr;
+ double abs_time_diff = 999;
+ for (auto sts_hit_j : fStsHits[sts_unit_j]) {
+
+ if (is_mc_data) {
+ int32_t sts_address_j = sts_hit_j->GetAddress();
+ double sts_j_x = sts_hit_j->GetX();
+ double sts_j_y = sts_hit_j->GetY();
+
+ fH2D[Form("Sts%d:Sts%d_YY", sts_unit_i, sts_unit_j)]->Fill(sts_i_y, sts_j_y);
+ fH2D[Form("Sts%d:Sts%d_XX", sts_unit_i, sts_unit_j)]->Fill(sts_i_x, sts_j_x);
+ fH2D[Form("Sts0x%x:Sts0x%x_YY", sts_address_i, sts_address_j)]->Fill(sts_i_y, sts_j_y);
+ fH2D[Form("Sts0x%x:Sts0x%x_XX", sts_address_i, sts_address_j)]->Fill(sts_i_x, sts_j_x);
+ }
+ else {
+
+ double dt = std::abs(sts_hit_i->GetTime() - sts_hit_j->GetTime());
+ if (abs_time_diff > dt) {
+ abs_time_diff = dt;
+ closest_hit = sts_hit_j;
+ }
+ }
+ }
+ if (!is_mc_data && closest_hit != nullptr) {
+ int32_t sts_address_j = closest_hit->GetAddress();
+ double sts_j_x = closest_hit->GetX();
+ double sts_j_y = closest_hit->GetY();
+
+ fH2D[Form("Sts%d:Sts%d_YY", sts_unit_i, sts_unit_j)]->Fill(sts_i_y, sts_j_y);
+ fH2D[Form("Sts%d:Sts%d_XX", sts_unit_i, sts_unit_j)]->Fill(sts_i_x, sts_j_x);
+ fH2D[Form("Sts0x%x:Sts0x%x_YY", sts_address_i, sts_address_j)]->Fill(sts_i_y, sts_j_y);
+ fH2D[Form("Sts0x%x:Sts0x%x_XX", sts_address_i, sts_address_j)]->Fill(sts_i_x, sts_j_x);
+ }
+ }
+ }
+ }
+ fStsHits.clear();
+}
+
+void CbmStsCorrelation::ProcessEvent(CbmEvent* evt)
+{
+ if (fAnalysisCuts != nullptr && !fAnalysisCuts->CheckEvent(evt)) return;
+
+ int nb_sts_hits = evt->GetNofData(ECbmDataType::kStsHit);
+ LOG(debug) << Form("Processing event with %d StsHit", nb_sts_hits);
+ for (int hit_evt_idx = 0; hit_evt_idx < nb_sts_hits; hit_evt_idx++) {
+ int sts_hit_idx = evt->GetIndex(ECbmDataType::kStsHit, hit_evt_idx);
+ CbmStsHit* hit = (CbmStsHit*) fStsHitArray->UncheckedAt(sts_hit_idx);
+ ProcessHit(hit);
+ }
+}
+
+void CbmStsCorrelation::ProcessHit(CbmStsHit* hit)
+{
+ if (hit == nullptr) return;
+
+ if (fAnalysisCuts != nullptr && !fAnalysisCuts->CheckStsHit(hit, fStsCluArray)) return;
+
+ int32_t address = hit->GetAddress();
+ int32_t unit = CbmStsAddress::GetElementId(address, kStsUnit);
+
+ fStsHits[unit].push_back(hit);
+}
+
+void CbmStsCorrelation::Exec(Option_t*)
+{
+ // Check if CbmEvent
+ if (fCbmEvtArray != nullptr) {
+ uint32_t nb_events = fCbmEvtArray->GetEntriesFast();
+ LOG(info) << Form("Running event-like Entry: %d\t%d CbmEvents", entry_, nb_events);
+ for (uint32_t evt_idx = 0; evt_idx < nb_events; evt_idx++) {
+ CbmEvent* event = (CbmEvent*) fCbmEvtArray->UncheckedAt(evt_idx);
+ ProcessEvent(event);
+ BuildCorrelation();
+ }
+ }
+ else {
+
+ LOG(debug) << "Running time-like ...";
+ uint32_t n_of_hits = fStsHitArray->GetEntriesFast();
+ for (uint32_t hit_idx = 0; hit_idx < n_of_hits; hit_idx++) {
+ CbmStsHit* sts_hit = (CbmStsHit*) fStsHitArray->UncheckedAt(hit_idx);
+ ProcessHit(sts_hit);
+ }
+ BuildCorrelation();
+ }
+ entry_++;
+}
+
+InitStatus CbmStsCorrelation::Init()
+{
+ LOG(debug) << "Init CbmStsCorrelation ...";
+
+ FairRootManager* ioman = FairRootManager::Instance();
+ if (ioman != nullptr) {
+ fCbmEvtArray = (TClonesArray*) ioman->GetObject("CbmEvent");
+ fStsHitArray = (TClonesArray*) ioman->GetObject("StsHit");
+ fStsCluArray = (TClonesArray*) ioman->GetObject("StsCluster");
+ }
+
+ LoadSetup();
+
+ BookHistograms();
+
+ return kSUCCESS;
+}
+
+
+void CbmStsCorrelation::Finish()
+{
+ SaveToFile();
+ fH1D.clear();
+ fH2D.clear();
+}
diff --git a/analysis/detectors/sts/CbmStsCorrelation.h b/analysis/detectors/sts/CbmStsCorrelation.h
new file mode 100644
index 0000000000000000000000000000000000000000..34cbf428cb28ac1a580bd36af7a5633a16da6d8e
--- /dev/null
+++ b/analysis/detectors/sts/CbmStsCorrelation.h
@@ -0,0 +1,65 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+#ifndef CBMSTSCORRELATION_H
+#define CBMSTSCORRELATION_H
+
+#include "CbmEvent.h"
+#include "CbmStsAnaBase.h"
+#include "CbmStsCluster.h"
+#include "CbmStsHit.h"
+#include "CbmStsUtils.h"
+
+#include <FairTask.h>
+
+#include <TClonesArray.h>
+
+#include <cstring>
+#include <map>
+#include <unordered_map>
+#include <unordered_set>
+
+/** \class CbmStsCorrelation
+ * \brief Task for correlation analysis of STS hits.
+ *
+ * This class inherits from FairTask and CbmStsAnaBase.
+ * It provides functionality for analyzing correlations
+ * among STS hits.
+ */
+class CbmStsCorrelation : public FairTask, public CbmStsAnaBase {
+ public:
+ CbmStsCorrelation() = default;
+ ~CbmStsCorrelation() = default;
+
+ InitStatus Init();
+ void Exec(Option_t*);
+ void Finish();
+
+ private:
+ std::map<int32_t, std::vector<CbmStsHit*>> fStsHits;
+
+ TClonesArray* fCbmEvtArray;
+ TClonesArray* fStsHitArray;
+ TClonesArray* fStsCluArray;
+
+ void BookHistograms();
+
+ /** \brief Process an Cbm events
+ * It filters event based on the provided CbmCutMap
+ */
+ void ProcessEvent(CbmEvent*);
+
+ /** \brief Process an STS hit
+ * It filters hits based on the provided CbmCutMap
+ */
+ void ProcessHit(CbmStsHit*);
+
+ /** \brief Build the correlation
+ * Using filtered hit build the combinatorial among different station sensors
+ */
+ void BuildCorrelation();
+
+ ClassDef(CbmStsCorrelation, 1);
+};
+#endif
diff --git a/analysis/detectors/sts/CbmStsEfficiency.cxx b/analysis/detectors/sts/CbmStsEfficiency.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..7823ea455b0ebe49a9272ad3f84ade8459a37cb5
--- /dev/null
+++ b/analysis/detectors/sts/CbmStsEfficiency.cxx
@@ -0,0 +1,583 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+#include "CbmStsEfficiency.h"
+
+#include "CbmStsAddress.h"
+#include "CbmStsUtils.h"
+
+CbmStsEfficiency::CbmStsEfficiency(uint32_t layer_idx) : fTestLayer(layer_idx)
+{
+ LOG(debug) << "Creating an instance of CbmStsEfficiency ...";
+}
+CbmStsEfficiency::CbmStsEfficiency(uint32_t test_layer, double max_dis_trg_vtx, double max_dca_trk_vtx,
+ double default_residual)
+ : fTestLayer(test_layer)
+ , fMaxDisTrgVtx(max_dis_trg_vtx)
+ , fMaxDCATrkVtx(max_dca_trk_vtx)
+ , fDefaultResidual(default_residual)
+{
+ LOG(debug) << "Creating an instance of CbmStsEfficiency ...";
+}
+
+void CbmStsEfficiency::BookHistograms()
+{
+ // ---- Binning configuration ----
+ // Vertex
+ double vertex_x_min = -15;
+ double vertex_x_max = +15;
+ double vertex_y_min = -15;
+ double vertex_y_max = +15;
+ double vertex_z_min = -15;
+ double vertex_z_max = +15;
+ double vertex_dx = 0.1;
+ double vertex_dy = 0.1;
+ double vertex_dz = 0.1;
+
+ auto x_vtx_binning =
+ cbm_sts_utils::HBinning{uint32_t((vertex_x_max - vertex_x_min) / vertex_dx), vertex_x_min, vertex_x_max};
+ auto y_vtx_binning =
+ cbm_sts_utils::HBinning{uint32_t((vertex_y_max - vertex_y_min) / vertex_dy), vertex_y_min, vertex_y_max};
+ auto z_vtx_binning =
+ cbm_sts_utils::HBinning{uint32_t((vertex_z_max - vertex_z_min) / vertex_dz), vertex_z_min, vertex_z_max};
+
+ // DCA
+ double dca_dx = .001; // 10 um
+ double dca_min = -2.5 - 0.5 * dca_dx;
+ double dca_max = +2.5 + 0.5 * dca_dx;
+ auto r_dca_binning = cbm_sts_utils::HBinning{uint32_t((dca_max - dca_min) / dca_dx), dca_min, dca_max};
+
+ // Residuals track - hit
+ double res_dx = +0.001; // 10 um
+ double res_min = -2.50 - 0.5 * res_dx;
+ double res_max = +2.50 + 0.5 * res_dx;
+ auto r_sen_binning = cbm_sts_utils::HBinning{uint32_t((res_max - res_min) / res_dx), res_min, res_max};
+ // ---- --------------------- ----
+
+ std::string h_name;
+ for (auto& [element, geo] : fStsGeoInfo) {
+ std::string h_name_base = "";
+ if (element == fTestLayer) {
+ h_name_base = Form("Sts%d", element);
+ }
+ else if (element > 8) {
+ if (CbmStsAddress::GetElementId(element, kStsUnit) != fTestLayer) continue;
+ fDUTList.push_back(element);
+ h_name_base = Form("Sts0x%x", element);
+ }
+
+ if (!h_name_base.length()) continue;
+
+ auto x_sen_binning = cbm_sts_utils::HBinning{uint32_t((geo[1] - geo[0]) / cbm_sts_utils::kStsDx), geo[0], geo[1]};
+ auto y_sen_binning = cbm_sts_utils::HBinning{uint32_t((geo[3] - geo[2]) / cbm_sts_utils::kStsDy), geo[2], geo[3]};
+
+ // ---- efficiency ----
+ h_name = Form("%s_found", h_name_base.c_str());
+ fH2D[h_name] =
+ std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), x_sen_binning.n_of_bins, x_sen_binning.x_min,
+ x_sen_binning.x_max, y_sen_binning.n_of_bins, y_sen_binning.x_min, y_sen_binning.x_max);
+ fH2D[h_name]->GetXaxis()->SetTitle("StsHit^{rec}_{X} [cm]");
+ fH2D[h_name]->GetYaxis()->SetTitle("StsHit^{rec}_{Y} [cm]");
+
+ h_name = Form("%s_track", h_name_base.c_str());
+ fH2D[h_name] =
+ std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), x_sen_binning.n_of_bins, x_sen_binning.x_min,
+ x_sen_binning.x_max, y_sen_binning.n_of_bins, y_sen_binning.x_min, y_sen_binning.x_max);
+ fH2D[h_name]->GetXaxis()->SetTitle("StsHit^{ext}_{X} [cm]");
+ fH2D[h_name]->GetYaxis()->SetTitle("StsHit^{ext}_{Y} [cm]");
+
+ h_name = Form("%s_hit_map", h_name_base.c_str());
+ fH2D[h_name] =
+ std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), x_sen_binning.n_of_bins, x_sen_binning.x_min,
+ x_sen_binning.x_max, y_sen_binning.n_of_bins, y_sen_binning.x_min, y_sen_binning.x_max);
+ fH2D[h_name]->GetXaxis()->SetTitle("StsHit_{X} [cm]");
+ fH2D[h_name]->GetYaxis()->SetTitle("StsHit_{Y} [cm]");
+ // ---- ---------- ----
+
+ // ---- residuals ----
+ for (const char* di : {"x", "y"}) {
+ h_name = Form("%s_d%s_vs_x", h_name_base.c_str(), di);
+ fH2D[h_name] =
+ std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), r_sen_binning.n_of_bins, r_sen_binning.x_min,
+ r_sen_binning.x_max, x_sen_binning.n_of_bins, x_sen_binning.x_min, x_sen_binning.x_max);
+ fH2D[h_name]->GetYaxis()->SetTitle("X [cm]");
+ fH2D[h_name]->GetXaxis()->SetTitle(Form("\\rho_{%s} [cm]", di));
+
+ h_name = Form("%s_d%s_vs_y", h_name_base.c_str(), di);
+ fH2D[h_name] =
+ std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), r_sen_binning.n_of_bins, r_sen_binning.x_min,
+ r_sen_binning.x_max, y_sen_binning.n_of_bins, y_sen_binning.x_min, y_sen_binning.x_max);
+ fH2D[h_name]->GetYaxis()->SetTitle("Y [cm]");
+ fH2D[h_name]->GetXaxis()->SetTitle(Form("\\rho_{%s} [cm]", di));
+ }
+ // ---- --------- ----
+ } // end geo loop
+
+ // ---- 3D vertex ----
+ h_name = "vertex_y_vs_x";
+ fH2D[h_name] =
+ std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), x_vtx_binning.n_of_bins, x_vtx_binning.x_min,
+ x_vtx_binning.x_max, y_vtx_binning.n_of_bins, y_vtx_binning.x_min, y_vtx_binning.x_max);
+ fH2D[h_name]->GetXaxis()->SetTitle("Vertex_{X} [cm]");
+ fH2D[h_name]->GetYaxis()->SetTitle("Vertex_{Y} [cm]");
+
+ h_name = "vertex_x_vs_z";
+ fH2D[h_name] =
+ std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), z_vtx_binning.n_of_bins, z_vtx_binning.x_min,
+ z_vtx_binning.x_max, x_vtx_binning.n_of_bins, x_vtx_binning.x_min, x_vtx_binning.x_max);
+ fH2D[h_name]->GetXaxis()->SetTitle("Vertex_{Z} [cm]");
+ fH2D[h_name]->GetYaxis()->SetTitle("Vertex_{X} [cm]");
+
+ h_name = "vertex_y_vs_z";
+ fH2D[h_name] =
+ std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), z_vtx_binning.n_of_bins, z_vtx_binning.x_min,
+ z_vtx_binning.x_max, y_vtx_binning.n_of_bins, y_vtx_binning.x_min, y_vtx_binning.x_max);
+ fH2D[h_name]->GetXaxis()->SetTitle("Vertex_{Z} [cm]");
+ fH2D[h_name]->GetYaxis()->SetTitle("Vertex_{Y} [cm]");
+ // ---- ------------------ ----
+
+ // ---- DCA ----
+ for (const char* di : {"x", "y", "z"}) {
+ h_name = Form("dca_d%s_vs_x", di);
+ fH2D[h_name] =
+ std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), r_dca_binning.n_of_bins, r_dca_binning.x_min,
+ r_dca_binning.x_max, x_vtx_binning.n_of_bins, x_vtx_binning.x_min, x_vtx_binning.x_max);
+ fH2D[h_name]->GetYaxis()->SetTitle("X [cm]");
+ fH2D[h_name]->GetXaxis()->SetTitle(Form("DCA_{%s} [cm]", di));
+
+ h_name = Form("dca_d%s_vs_y", di);
+ fH2D[h_name] =
+ std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), r_dca_binning.n_of_bins, r_dca_binning.x_min,
+ r_dca_binning.x_max, y_vtx_binning.n_of_bins, y_vtx_binning.x_min, y_vtx_binning.x_max);
+ fH2D[h_name]->GetYaxis()->SetTitle("Y [cm]");
+ fH2D[h_name]->GetXaxis()->SetTitle(Form("DCA_{%s} [cm]", di));
+
+ h_name = Form("dca_d%s_vs_z", di);
+ fH2D[h_name] =
+ std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), r_dca_binning.n_of_bins, r_dca_binning.x_min,
+ r_dca_binning.x_max, z_vtx_binning.n_of_bins, z_vtx_binning.x_min, z_vtx_binning.x_max);
+ fH2D[h_name]->GetYaxis()->SetTitle(Form("Z [cm]"));
+ fH2D[h_name]->GetXaxis()->SetTitle(Form("DCA_{%s} [cm]", di));
+ }
+ // ---- ------------------ ----
+
+ // ---- Track multiplicity ----
+ for (const char* mod : {":all", ":sel"}) {
+ h_name = Form("ca_track_multiplicity%s", mod);
+ fH1D[h_name] = std::make_unique<TH1D>(h_name.c_str(), h_name.c_str(), 50, 0, 50); // HARDCODE
+ fH1D[h_name]->GetXaxis()->SetTitle("N_{CATrack}");
+ fH1D[h_name]->GetYaxis()->SetTitle("Entries");
+
+ h_name = Form("ca_track_chi2%s", mod);
+ fH1D[h_name] = std::make_unique<TH1D>(h_name.c_str(), h_name.c_str(), 10000, 0, 10); // HARDCODE
+ fH1D[h_name]->GetXaxis()->SetTitle("Chi2");
+ fH1D[h_name]->GetYaxis()->SetTitle("Entries");
+ }
+ // ---- ------------------ ----
+}
+
+void CbmStsEfficiency::CheckEfficiency()
+{
+ LOG(debug) << "Checking efficiency of STS station " << fTestLayer << " using " << fGlbTrks.size()
+ << " CbmGlobalTracks";
+ fVertexFinder->SetTracks(fGlbTrks);
+ const std::optional<CbmVertex> vertex = fVertexFinder->FindVertex();
+
+ if (!vertex.has_value()) return;
+
+ if (std::abs(vertex->GetZ()) > fMaxDisTrgVtx) return;
+
+ fH2D["vertex_y_vs_x"]->Fill(vertex->GetX(), vertex->GetY());
+ fH2D["vertex_x_vs_z"]->Fill(vertex->GetZ(), vertex->GetX());
+ fH2D["vertex_y_vs_z"]->Fill(vertex->GetZ(), vertex->GetY());
+
+ TVector3 vtx(vertex->GetX(), vertex->GetY(), vertex->GetZ());
+ std::map<int32_t, std::vector<bool>> hit_used;
+ for (auto trk : fGlbTrks) {
+ double trk_tx = trk->GetParamFirst()->GetTx();
+ double trk_ty = trk->GetParamFirst()->GetTy();
+ double trk_x0 = trk->GetParamFirst()->GetX();
+ double trk_y0 = trk->GetParamFirst()->GetY();
+ double trk_z0 = trk->GetParamFirst()->GetZ();
+
+ // Check if tack comes from vertex
+ TVector3 p(trk_x0, trk_y0, trk_z0);
+ TVector3 e(trk_tx, trk_ty, 1);
+ TVector3 trk_to_vtx = ((vtx - p).Cross(e)) * (1. / e.Mag());
+
+ fH2D["dca_dx_vs_x"]->Fill(trk_to_vtx.Px(), vtx.Px());
+ fH2D["dca_dx_vs_y"]->Fill(trk_to_vtx.Px(), vtx.Py());
+ fH2D["dca_dx_vs_z"]->Fill(trk_to_vtx.Px(), vtx.Pz());
+
+ fH2D["dca_dy_vs_x"]->Fill(trk_to_vtx.Py(), vtx.Px());
+ fH2D["dca_dy_vs_y"]->Fill(trk_to_vtx.Py(), vtx.Py());
+ fH2D["dca_dy_vs_z"]->Fill(trk_to_vtx.Py(), vtx.Pz());
+
+ fH2D["dca_dz_vs_x"]->Fill(trk_to_vtx.Pz(), vtx.Px());
+ fH2D["dca_dz_vs_y"]->Fill(trk_to_vtx.Pz(), vtx.Py());
+ fH2D["dca_dz_vs_z"]->Fill(trk_to_vtx.Pz(), vtx.Pz());
+
+ if (trk_to_vtx.Mag() > fMaxDCATrkVtx) continue;
+
+ for (int32_t& sensor : fDUTList) {
+ auto hits = fStsHits[sensor];
+ size_t n_of_hits = hits.size();
+
+ // Extrapolate to the current sensor z-plane
+ double sensor_z = 0.5 * (fStsGeoInfo[sensor][4] + fStsGeoInfo[sensor][5]);
+ double predicted_x = trk_x0 + trk_tx * (sensor_z - trk_z0);
+ double predicted_y = trk_y0 + trk_ty * (sensor_z - trk_z0);
+
+ // Check if the sensor contain the space point
+ if (predicted_x < fStsGeoInfo[sensor][0] || predicted_x > fStsGeoInfo[sensor][1]
+ || predicted_y < fStsGeoInfo[sensor][2] || predicted_y > fStsGeoInfo[sensor][3]) {
+ continue;
+ }
+
+ int unit = CbmStsAddress::GetElementId(sensor, kStsUnit);
+ fH2D[Form("Sts%d_track", unit)]->Fill(predicted_x, predicted_y);
+ fH2D[Form("Sts0x%x_track", sensor)]->Fill(predicted_x, predicted_y);
+
+ if (n_of_hits) {
+ /* Block to search in available hits */
+ double closest_id = 0;
+ double dx = hits[0]->GetX() - predicted_x - fResiduals[sensor].mean.Px();
+ double dy = hits[0]->GetY() - predicted_y - fResiduals[sensor].mean.Py();
+ double closest_rho = dx * dx + dy * dy;
+
+ for (size_t idx = 1; idx < n_of_hits; idx++) {
+ dx = hits[idx]->GetX() - predicted_x - fResiduals[sensor].mean.Px();
+ dy = hits[idx]->GetY() - predicted_y - fResiduals[sensor].mean.Py();
+ double rho = dx * dx + dy * dy;
+
+ if (rho < closest_rho) {
+ closest_id = idx;
+ closest_rho = rho;
+ }
+ }
+
+ dx = hits[closest_id]->GetX() - predicted_x - fResiduals[sensor].mean.Px();
+ dy = hits[closest_id]->GetY() - predicted_y - fResiduals[sensor].mean.Py();
+
+ fH2D[Form("Sts%d_dx_vs_x", unit)]->Fill(dx, hits[closest_id]->GetX());
+ fH2D[Form("Sts%d_dy_vs_x", unit)]->Fill(dy, hits[closest_id]->GetX());
+ fH2D[Form("Sts%d_dx_vs_y", unit)]->Fill(dx, hits[closest_id]->GetY());
+ fH2D[Form("Sts%d_dy_vs_y", unit)]->Fill(dy, hits[closest_id]->GetY());
+
+ fH2D[Form("Sts0x%x_dx_vs_x", sensor)]->Fill(dx, hits[closest_id]->GetX());
+ fH2D[Form("Sts0x%x_dy_vs_x", sensor)]->Fill(dy, hits[closest_id]->GetX());
+ fH2D[Form("Sts0x%x_dx_vs_y", sensor)]->Fill(dx, hits[closest_id]->GetY());
+ fH2D[Form("Sts0x%x_dy_vs_y", sensor)]->Fill(dy, hits[closest_id]->GetY());
+
+ double sensor_resolution =
+ fResiduals[sensor].sigma.Mag() != 0 ? fResiduals[sensor].sigma.Mag() : fDefaultResidual;
+ if (closest_rho <= 3.0 * sensor_resolution) {
+ fH2D[Form("Sts%d_found", unit)]->Fill(predicted_x, predicted_y);
+ fH2D[Form("Sts0x%x_found", sensor)]->Fill(predicted_x, predicted_y);
+ }
+ }
+
+
+ /* Block to search in avaliable hit */
+
+ } // end sensor under test lopp
+ } // end track loop
+}
+
+TH2D* FindInactiveArea(TH2D* h, int dead_size_x = 10, int dead_size_y = 10)
+{
+ TH2D* dead = (TH2D*) h->Clone();
+ dead->Reset();
+
+ for (int bin_y = 1; bin_y < dead->GetNbinsY(); bin_y++) {
+ for (int bin_x = dead_size_x; bin_x < dead->GetNbinsX() - dead_size_x; bin_x++) {
+
+ bool all_dead = true;
+ for (int test_bin = -dead_size_x; test_bin <= +dead_size_x; test_bin++) {
+ if (h->GetBinContent(bin_x + test_bin, bin_y)) {
+ all_dead = false;
+ break;
+ }
+ }
+ if (all_dead) {
+ for (int test_bin = -dead_size_x; test_bin <= +dead_size_x; test_bin++)
+ dead->SetBinContent(bin_x + test_bin, bin_y, 1);
+ }
+ }
+ }
+
+ for (int bin_x = 1; bin_x < dead->GetNbinsX(); bin_x++) {
+ for (int bin_y = dead_size_y; bin_y < dead->GetNbinsY() - dead_size_y; bin_y++) {
+
+ bool all_dead = true;
+ for (int test_bin = -dead_size_y; test_bin <= +dead_size_y; test_bin++) {
+ if (h->GetBinContent(bin_x, bin_y + test_bin)) {
+ all_dead = false;
+ break;
+ }
+ }
+ if (all_dead) {
+ for (int test_bin = -dead_size_y; test_bin <= +dead_size_y; test_bin++)
+ dead->SetBinContent(bin_x, bin_y + test_bin, 1);
+ }
+ }
+ }
+
+ return dead;
+}
+
+void CbmStsEfficiency::Efficiency(uint32_t sensor)
+{
+ std::string h_name_base = "";
+ if (sensor == fTestLayer) {
+ h_name_base = Form("Sts%d", sensor);
+ }
+ else if (sensor > 8) {
+ if (CbmStsAddress::GetElementId(sensor, kStsUnit) != fTestLayer) return;
+ h_name_base = Form("Sts0x%x", sensor);
+ }
+ auto found = fH2D[Form("%s_found", h_name_base.c_str())].get();
+ auto track = fH2D[Form("%s_track", h_name_base.c_str())].get();
+
+ if (found == nullptr || track == nullptr) return;
+
+ TH2D* inactive_area = FindInactiveArea(fH2D[Form("%s_hit_map", h_name_base.c_str())].get(), 50, 50);
+
+ // Manual integration
+ int n_of_found = 0;
+ int n_of_track = 0;
+ std::vector<double> samples;
+
+ for (int bin_y = 50; bin_y < inactive_area->GetNbinsY() - 50; bin_y++) {
+ for (int bin_x = 50; bin_x < inactive_area->GetNbinsX() - 50; bin_x++) {
+ if (!inactive_area->GetBinContent(bin_x, bin_y)) {
+ n_of_found += found->GetBinContent(bin_x, bin_y);
+ n_of_track += track->GetBinContent(bin_x, bin_y);
+ }
+ }
+ }
+ LOG(debug) << n_of_found << "\t" << n_of_track;
+ double efficiency = n_of_track != 0 ? double(n_of_found) / n_of_track : -1;
+ if (efficiency != -1) {
+ samples.push_back(efficiency);
+ }
+
+ double mean = 0;
+ double sigma = 0;
+ int n = samples.size();
+ if (n == 0) {
+ return;
+ }
+
+ // Take efficiency as mean value
+ for (double& val : samples) {
+ mean += val;
+ }
+ mean = double(mean / n);
+
+ // Take error as RMS of the values
+ for (double& val : samples) {
+ sigma += (val - mean) * (val - mean);
+ }
+ sigma = sqrt(double(sigma / n));
+
+ std::cout << Form("0x%x:\t%0.4f\t%0.4f\n", sensor, mean, sigma);
+}
+
+void CbmStsEfficiency::FinishTask()
+{
+ for (auto& [element, geo] : fStsGeoInfo) {
+ Efficiency(element);
+ }
+
+ SaveToFile();
+}
+
+
+void CbmStsEfficiency::ProcessEvent(CbmEvent* evt)
+{
+ int nb_sts_hits = evt->GetNofData(ECbmDataType::kStsHit);
+
+ for (int hit_evt_idx = 0; hit_evt_idx < nb_sts_hits; hit_evt_idx++) {
+ int sts_hit_idx = evt->GetIndex(ECbmDataType::kStsHit, hit_evt_idx);
+ CbmStsHit* hit = (CbmStsHit*) fStsHitArray->UncheckedAt(sts_hit_idx);
+ ProcessHit(hit);
+ }
+
+ int nb_of_glob_trk = evt->GetNofData(ECbmDataType::kGlobalTrack);
+ for (int evt_glob_trk_idx = 0; evt_glob_trk_idx < nb_of_glob_trk; evt_glob_trk_idx++) {
+ int glob_trk_idx = evt->GetIndex(ECbmDataType::kGlobalTrack, evt_glob_trk_idx);
+ CbmGlobalTrack* glob_trk = (CbmGlobalTrack*) fGlbTrkArray->UncheckedAt(glob_trk_idx);
+ ProcessGlobalTrack(glob_trk);
+ }
+
+ fH1D["ca_track_multiplicity:all"]->Fill(nb_of_glob_trk);
+ fH1D["ca_track_multiplicity:sel"]->Fill(fGlbTrks.size());
+
+ if (fGlbTrks.size()) CheckEfficiency();
+
+ fGlbTrks.clear();
+ fStsHits.clear();
+}
+
+void CbmStsEfficiency::ProcessGlobalTrack(CbmGlobalTrack* trk)
+{
+ auto sts_trk_idx = trk->GetStsTrackIndex();
+ auto rich_trk_idx = trk->GetRichRingIndex();
+ auto much_trk_idx = trk->GetMuchTrackIndex();
+ auto trd_trk_idx = trk->GetTrdTrackIndex();
+ auto tof_trk_idx = trk->GetTofTrackIndex();
+ double trk_p_val = TMath::Prob(trk->GetChi2(), trk->GetNDF());
+ fH1D["ca_track_chi2:all"]->Fill(trk->GetChi2());
+
+ // Apply GlobalTracks cuts
+ int32_t mvd_trk_size = sts_trk_idx != -1 && fStsTrkArray != nullptr
+ ? ((CbmStsTrack*) fStsTrkArray->UncheckedAt(sts_trk_idx))->GetNofMvdHits()
+ : 0;
+
+ int32_t sts_trk_size = sts_trk_idx != -1 && fStsTrkArray != nullptr
+ ? ((CbmStsTrack*) fStsTrkArray->UncheckedAt(sts_trk_idx))->GetNofStsHits()
+ : 0;
+
+ int32_t rich_trk_size = rich_trk_idx != -1 && fRchTrkArray != nullptr
+ ? ((CbmTrack*) fRchTrkArray->UncheckedAt(rich_trk_idx))->GetNofHits()
+ : 0;
+
+ int32_t much_trk_size = much_trk_idx != -1 && fMchTrkArray != nullptr
+ ? ((CbmTrack*) fMchTrkArray->UncheckedAt(much_trk_idx))->GetNofHits()
+ : 0;
+
+ int32_t trd_trk_size = trd_trk_idx != -1 && fTrdTrkArray != nullptr
+ ? ((CbmTrack*) fTrdTrkArray->UncheckedAt(trd_trk_idx))->GetNofHits()
+ : 0;
+
+ int32_t tof_trk_size = tof_trk_idx != -1 && fTofTrkArray != nullptr
+ ? ((CbmTrack*) fTofTrkArray->UncheckedAt(tof_trk_idx))->GetNofHits()
+ : 0;
+
+ if (fAnalysisCuts != nullptr
+ && (!fAnalysisCuts->Check(CbmCutId::kGlobalTrackMvdSize, mvd_trk_size)
+ || !fAnalysisCuts->Check(CbmCutId::kGlobalTrackStsSize, sts_trk_size)
+ || !fAnalysisCuts->Check(CbmCutId::kGlobalTrackRichSize, rich_trk_size)
+ || !fAnalysisCuts->Check(CbmCutId::kGlobalTrackMuchSize, much_trk_size)
+ || !fAnalysisCuts->Check(CbmCutId::kGlobalTrackTrdSize, trd_trk_size)
+ || !fAnalysisCuts->Check(CbmCutId::kGlobalTrackTofSize, tof_trk_size)
+ || !fAnalysisCuts->Check(CbmCutId::kGlobalTrackChi2, trk->GetChi2())
+ || !fAnalysisCuts->Check(CbmCutId::kGlobalTrackPval, trk_p_val))) {
+ return;
+ }
+
+ fH1D["ca_track_chi2:sel"]->Fill(trk->GetChi2());
+ fGlbTrks.push_back(trk);
+}
+
+void CbmStsEfficiency::ProcessHit(CbmStsHit* hit)
+{
+ if (hit == nullptr) return;
+ uint32_t address = hit->GetAddress();
+ uint32_t unit = CbmStsAddress::GetElementId(address, kStsUnit);
+
+ if (unit != fTestLayer) return;
+
+ if (fAnalysisCuts != nullptr && !fAnalysisCuts->CheckStsHit(hit, fStsCluArray)) return;
+
+ fStsHits[address].push_back(hit);
+ fH2D[Form("Sts%u_hit_map", unit)]->Fill(hit->GetX(), hit->GetY());
+ fH2D[Form("Sts0x%x_hit_map", address)]->Fill(hit->GetX(), hit->GetY());
+}
+
+void CbmStsEfficiency::Exec(Option_t*)
+{
+ // Check if CbmEvent
+ if (fCbmEvtArray != nullptr) {
+ LOG(info) << "Running CbmStsEfficiency - Event-like - " << entry_;
+
+ uint32_t nb_events = fCbmEvtArray->GetEntriesFast();
+ for (uint32_t evt_idx = 0; evt_idx < nb_events; evt_idx++) {
+ CbmEvent* event = (CbmEvent*) fCbmEvtArray->UncheckedAt(evt_idx);
+ if (fAnalysisCuts != nullptr && !fAnalysisCuts->CheckEvent(event)) continue;
+ LOG(debug) << Form("Process CbmEvent %d / %d", evt_idx, nb_events);
+ ProcessEvent(event);
+ }
+ }
+ else {
+
+ LOG(info) << "Running CbmStsEfficiency - TS-like - ";
+ uint32_t n_of_hits = fStsHitArray->GetEntriesFast();
+ for (uint32_t hit_idx = 0; hit_idx < n_of_hits; hit_idx++) {
+ CbmStsHit* sts_hit = (CbmStsHit*) fStsHitArray->UncheckedAt(hit_idx);
+ ProcessHit(sts_hit);
+ }
+ uint32_t nb_of_glob_trk = fGlbTrkArray->GetEntriesFast();
+ LOG(debug) << Form("Number of GlobalTracks: %u", nb_of_glob_trk);
+ for (uint32_t glob_trk_idx = 0; glob_trk_idx < nb_of_glob_trk; glob_trk_idx++) {
+ CbmGlobalTrack* glob_trk = (CbmGlobalTrack*) fGlbTrkArray->UncheckedAt(glob_trk_idx);
+ ProcessGlobalTrack(glob_trk);
+ }
+
+ if (fGlbTrks.size()) CheckEfficiency();
+
+ fGlbTrks.clear();
+ fStsHits.clear();
+ }
+ entry_++;
+}
+
+InitStatus CbmStsEfficiency::Init()
+{
+ LOG(debug) << "Init CbmStsEfficiency ...";
+
+ FairRootManager* ioman = FairRootManager::Instance();
+ if (ioman != nullptr) {
+ fCbmEvtArray = (TClonesArray*) ioman->GetObject("CbmEvent");
+
+ fGlbTrkArray = (TClonesArray*) ioman->GetObject("GlobalTrack");
+
+ fStsTrkArray = (TClonesArray*) ioman->GetObject("StsTrack");
+ fRchTrkArray = (TClonesArray*) ioman->GetObject("RichTrack");
+ fMchTrkArray = (TClonesArray*) ioman->GetObject("MuchTrack");
+ fTrdTrkArray = (TClonesArray*) ioman->GetObject("TrdTrack");
+ fTofTrkArray = (TClonesArray*) ioman->GetObject("TofTrack");
+
+ fStsHitArray = (TClonesArray*) ioman->GetObject("StsHit");
+
+ fStsCluArray = (TClonesArray*) ioman->GetObject("StsCluster");
+ }
+
+ LoadSetup();
+
+ BookHistograms();
+
+ std::stringstream ss;
+ for (auto s : fDUTList) {
+ ss << Form("\t\t0x%x - U%d L%d M%d\n", s, CbmStsAddress::GetElementId(s, kStsUnit),
+ CbmStsAddress::GetElementId(s, kStsLadder), CbmStsAddress::GetElementId(s, kStsModule));
+ }
+
+ LOG(info) << "Max vtx target dz: " << fMaxDisTrgVtx;
+ LOG(info) << "Max trk vtx DCA: " << fMaxDCATrkVtx;
+ LOG(info) << "Default residual: " << fDefaultResidual;
+ LOG(info) << "Test layer: " << fTestLayer;
+ LOG(info) << "Layer modules:\n" << ss.str();
+
+ return kSUCCESS;
+}
+
+void CbmStsEfficiency::SetSensorResidual(int32_t address, TVector3 m, TVector3 s)
+{
+ fResiduals[address] = Residual(m, s);
+}
+
+void CbmStsEfficiency::SetResidual(std::string file_name)
+{
+ std::ifstream in(file_name);
+ int32_t sensor;
+ double mean_x, mean_y, sigma_x, sigma_y;
+ while (in >> std::hex >> sensor >> std::dec >> mean_x >> mean_y >> sigma_x >> sigma_y) {
+ SetSensorResidual(sensor, TVector3(mean_x, mean_y, 0), TVector3(sigma_x, sigma_y, 0));
+ }
+}
+
+
+void CbmStsEfficiency::SetVertexFinder(CbmDcaVertexFinder* vtx_finder) { fVertexFinder = vtx_finder; }
diff --git a/analysis/detectors/sts/CbmStsEfficiency.h b/analysis/detectors/sts/CbmStsEfficiency.h
new file mode 100644
index 0000000000000000000000000000000000000000..29469fb0627748841808130e911e3b62ab82f7e8
--- /dev/null
+++ b/analysis/detectors/sts/CbmStsEfficiency.h
@@ -0,0 +1,115 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+#ifndef CBMSTSEFFICIENCY_H
+#define CBMSTSEFFICIENCY_H
+
+#include "CbmDcaVertexFinder.h"
+#include "CbmEvent.h"
+#include "CbmGlobalTrack.h"
+#include "CbmStsAnaBase.h"
+#include "CbmStsCluster.h"
+#include "CbmStsHit.h"
+#include "CbmStsTrack.h"
+#include "CbmStsUtils.h"
+#include "CbmTrack.h"
+
+#include <FairTask.h>
+
+#include <TClonesArray.h>
+#include <TVector3.h>
+
+#include <cstring>
+#include <map>
+#include <unordered_map>
+#include <unordered_set>
+
+struct Residual {
+ TVector3 mean;
+ TVector3 sigma;
+ Residual() {}
+ Residual(TVector3 m, TVector3 s)
+ {
+ mean = m;
+ sigma = s;
+ }
+};
+
+/** \class CbmStsEfficiency
+ * \brief Task for Hit Reconstruction Efficiency (HRE) analysis of STS hits.
+ *
+ * This class inherits from FairTask and CbmStsAnaBase.
+ * It provides functionality for analyzing the CA Tracking based HRE
+ * of STS hits.
+ */
+class CbmStsEfficiency : public FairTask, public CbmStsAnaBase {
+ public:
+ CbmStsEfficiency() = default;
+ ~CbmStsEfficiency() = default;
+
+ /** \brief Parameterized constructor */
+ CbmStsEfficiency(uint32_t);
+
+ /** \brief Parameterized constructor */
+ CbmStsEfficiency(uint32_t, double, double, double);
+
+ InitStatus Init();
+ void Exec(Option_t*);
+ void FinishTask();
+
+ void SetSensorResidual(int32_t, TVector3, TVector3);
+ void SetResidual(std::string);
+
+ void SetVertexFinder(CbmDcaVertexFinder*);
+
+ bool fDrawOpt{true};
+
+ private:
+ uint32_t fTestLayer{0};
+ std::vector<int32_t> fDUTList;
+ double fMaxDisTrgVtx{0.25};
+ double fMaxDCATrkVtx{0.25};
+ double fDefaultResidual{0.12};
+ std::map<int32_t, Residual> fResiduals;
+
+
+ CbmDcaVertexFinder* fVertexFinder;
+
+ std::vector<CbmGlobalTrack*> fGlbTrks;
+ std::map<int32_t, std::vector<CbmStsHit*>> fStsHits;
+
+ TClonesArray* fCbmEvtArray{nullptr};
+ TClonesArray* fGlbTrkArray{nullptr};
+ TClonesArray* fStsTrkArray{nullptr};
+ TClonesArray* fRchTrkArray{nullptr};
+ TClonesArray* fMchTrkArray{nullptr};
+ TClonesArray* fTrdTrkArray{nullptr};
+ TClonesArray* fTofTrkArray{nullptr};
+ TClonesArray* fStsHitArray{nullptr};
+ TClonesArray* fStsCluArray{nullptr};
+
+ void BookHistograms();
+
+ void Efficiency(uint32_t);
+
+ /** \brief Process an Cbm events
+ * It filters event based on the provided CbmCutMap
+ */
+ void ProcessEvent(CbmEvent*);
+
+ /** \brief Process an Global tracks
+ * It filters the tracks based on the provided CbmCutMap
+ */
+ void ProcessGlobalTrack(CbmGlobalTrack*);
+
+ /** \brief Process an STS hit
+ * It filters hits based on the provided CbmCutMap
+ */
+ void ProcessHit(CbmStsHit*);
+
+ void CheckEfficiency();
+
+ ClassDef(CbmStsEfficiency, 1);
+};
+#endif
diff --git a/analysis/detectors/sts/CbmStsHitAna.cxx b/analysis/detectors/sts/CbmStsHitAna.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..15785ec3d65274833296df70c9d1ff68a5b8b6e0
--- /dev/null
+++ b/analysis/detectors/sts/CbmStsHitAna.cxx
@@ -0,0 +1,241 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+#include "CbmStsHitAna.h"
+
+#include <CbmStsTrack.h>
+
+#include <ctime>
+#include <iostream>
+#include <typeinfo>
+
+CbmStsHitAna::CbmStsHitAna() : fHitModifier({":all"}) {}
+
+CbmStsHitAna::CbmStsHitAna(std::string cal_par_file) : fHitModifier({":all"}), fCalibrationFile(cal_par_file) {}
+
+InitStatus CbmStsHitAna::Init()
+{
+ if (fModuleParSet == nullptr) {
+ LOG(debug) << "Loading charge calibration ...";
+ fModuleParSet = std::make_unique<CbmStsParSetModule>("CbmParSetModule", "STS parameters", "Default");
+ fModuleParSet->LoadParASCII(fCalibrationFile);
+ }
+
+ FairRootManager* ioman = FairRootManager::Instance();
+ if (ioman == nullptr) return kERROR;
+
+ fCbmEvtArray = (TClonesArray*) ioman->GetObject("CbmEvent");
+
+ fGlbTrkArray = (TClonesArray*) ioman->GetObject("GlobalTrack");
+ fStsTrkArray = (TClonesArray*) ioman->GetObject("StsTrack");
+ fRchTrkArray = (TClonesArray*) ioman->GetObject("RichTrack");
+ fMchTrkArray = (TClonesArray*) ioman->GetObject("MuchTrack");
+ fTrdTrkArray = (TClonesArray*) ioman->GetObject("TrdTrack");
+ fTofTrkArray = (TClonesArray*) ioman->GetObject("TofTrack");
+
+ fStsHitArray = (TClonesArray*) ioman->GetObject("StsHit");
+ fStsCluArray = (TClonesArray*) ioman->GetObject("StsCluster");
+
+ if (fGlbTrkArray) fHitModifier.push_back(":trk");
+
+ LoadSetup();
+
+ return kSUCCESS;
+}
+
+void CbmStsHitAna::BookHistograms(int32_t address)
+{
+ LOG(debug) << Form("Booking histograms for module: 0x%x", address);
+
+ const CbmStsParModule& par_module = fModuleParSet->GetParModule(address);
+ const auto [n_side_binning, p_side_binning] = cbm_sts_utils::ChargeBinning(par_module, fMaxCluSize);
+
+ std::string h_name;
+
+ double x_min = fStsGeoInfo.count(address) ? fStsGeoInfo[address][0] : -15;
+ double x_max = fStsGeoInfo.count(address) ? fStsGeoInfo[address][1] : +15;
+ double y_min = fStsGeoInfo.count(address) ? fStsGeoInfo[address][2] : -15;
+ double y_max = fStsGeoInfo.count(address) ? fStsGeoInfo[address][3] : +15;
+
+ double t_min = -150 - 0.5 * cbm_sts_utils::kStsClock;
+ double t_max = +150 + 0.5 * cbm_sts_utils::kStsClock;
+
+ auto x_binning = cbm_sts_utils::HBinning{uint32_t((x_max - x_min) / cbm_sts_utils::kStsDx), x_min, x_max};
+ auto y_binning = cbm_sts_utils::HBinning{uint32_t((y_max - y_min) / cbm_sts_utils::kStsDy), y_min, y_max};
+ auto t_binning = cbm_sts_utils::HBinning{uint32_t((t_max - t_min) / cbm_sts_utils::kStsClock), t_min, t_max};
+
+ for (const char* mod : fHitModifier) {
+ h_name = Form("0x%x_Q_asymmetry%s", address, mod);
+ fH1D[h_name] = std::make_unique<TH1D>(h_name.c_str(), h_name.c_str(), 1000, -1.01, +1.01);
+
+ h_name = Form("0x%x_Qp_vs_Qn%s", address, mod);
+ fH2D[h_name] = std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), n_side_binning.n_of_bins,
+ n_side_binning.x_min, n_side_binning.x_max, p_side_binning.n_of_bins,
+ p_side_binning.x_min, p_side_binning.x_max);
+
+ h_name = Form("0x%x_Qp_vs_size%s", address, mod);
+ fH2D[h_name] = std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), p_side_binning.n_of_bins,
+ p_side_binning.x_min, p_side_binning.x_max, fMaxCluSize, 1, fMaxCluSize + 1);
+
+ h_name = Form("0x%x_Qn_vs_size%s", address, mod);
+ fH2D[h_name] = std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), n_side_binning.n_of_bins,
+ n_side_binning.x_min, n_side_binning.x_max, fMaxCluSize, 1, fMaxCluSize + 1);
+
+ h_name = Form("0x%x_psize_vs_nsize%s", address, mod);
+ fH2D[h_name] = std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), fMaxCluSize, 1, fMaxCluSize + 1, fMaxCluSize,
+ 1, fMaxCluSize + 1);
+
+ h_name = Form("0x%x_y_vs_x%s", address, mod);
+ fH2D[h_name] = std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), x_binning.n_of_bins, x_binning.x_min,
+ x_binning.x_max, y_binning.n_of_bins, y_binning.x_min, y_binning.x_max);
+
+ h_name = Form("0x%x_cluster_dt%s", address, mod);
+ fH1D[h_name] =
+ std::make_unique<TH1D>(h_name.c_str(), h_name.c_str(), t_binning.n_of_bins, t_binning.x_min, t_binning.x_max);
+ }
+}
+
+void CbmStsHitAna::ProcessEvent(CbmEvent* evt)
+{
+ if (fAnalysisCuts != nullptr && !fAnalysisCuts->CheckEvent(evt)) return;
+
+ int nb_sts_hits = evt->GetNofData(ECbmDataType::kStsHit);
+ for (int hit_evt_idx = 0; hit_evt_idx < nb_sts_hits; hit_evt_idx++) {
+ int sts_hit_idx = evt->GetIndex(ECbmDataType::kStsHit, hit_evt_idx);
+ CbmStsHit* hit = (CbmStsHit*) fStsHitArray->UncheckedAt(sts_hit_idx);
+ ProcessHit(hit, false);
+ }
+
+ int nb_of_glob_trk = evt->GetNofData(ECbmDataType::kGlobalTrack);
+ for (int evt_glob_trk_idx = 0; evt_glob_trk_idx < nb_of_glob_trk; evt_glob_trk_idx++) {
+ int glob_trk_idx = evt->GetIndex(ECbmDataType::kGlobalTrack, evt_glob_trk_idx);
+ CbmGlobalTrack* glob_trk = (CbmGlobalTrack*) fGlbTrkArray->UncheckedAt(glob_trk_idx);
+ ProcessGlobalTrack(glob_trk);
+ }
+}
+
+void CbmStsHitAna::ProcessHit(CbmStsHit* hit, bool belong_to_trk)
+{
+ if (hit == nullptr) return;
+ int32_t address = hit->GetAddress();
+
+ if (fAnalysisCuts != nullptr && !fAnalysisCuts->CheckStsHit(hit, fStsCluArray)) {
+ return;
+ }
+
+ if (!fAddressBook.count(address)) {
+ fAddressBook.insert(address);
+ BookHistograms(address);
+ }
+
+ double q_n = cbm_sts_utils::GetHitChargeF(hit, fStsCluArray);
+ double q_p = cbm_sts_utils::GetHitChargeB(hit, fStsCluArray);
+
+ int32_t size_n = cbm_sts_utils::GetHitCluSizeF(hit, fStsCluArray);
+ int32_t size_p = cbm_sts_utils::GetHitCluSizeB(hit, fStsCluArray);
+
+ const char* mod = fHitModifier[belong_to_trk];
+
+ fH1D[Form("0x%x_Q_asymmetry%s", address, mod)]->Fill(cbm_sts_utils::GetHitChargeAsy(hit, fStsCluArray));
+ fH2D[Form("0x%x_Qp_vs_Qn%s", address, mod)]->Fill(q_n, q_p);
+ fH2D[Form("0x%x_Qp_vs_size%s", address, mod)]->Fill(q_p, size_p);
+ fH2D[Form("0x%x_Qn_vs_size%s", address, mod)]->Fill(q_n, size_n);
+ fH2D[Form("0x%x_psize_vs_nsize%s", address, mod)]->Fill(size_p, size_n);
+ fH2D[Form("0x%x_y_vs_x%s", address, mod)]->Fill(hit->GetX(), hit->GetY());
+ fH1D[Form("0x%x_cluster_dt%s", address, mod)]->Fill(cbm_sts_utils::GetHitTimeF(hit, fStsCluArray)
+ - cbm_sts_utils::GetHitTimeB(hit, fStsCluArray));
+}
+
+void CbmStsHitAna::ProcessGlobalTrack(CbmGlobalTrack* trk)
+{
+ auto sts_trk_idx = trk->GetStsTrackIndex();
+ auto rich_trk_idx = trk->GetRichRingIndex();
+ auto much_trk_idx = trk->GetMuchTrackIndex();
+ auto trd_trk_idx = trk->GetTrdTrackIndex();
+ auto tof_trk_idx = trk->GetTofTrackIndex();
+ float trk_p_val = TMath::Prob(trk->GetChi2(), trk->GetNDF());
+
+ // Apply GlobalTracks cuts
+ int32_t mvd_trk_size = sts_trk_idx != -1 && fStsTrkArray != nullptr
+ ? ((CbmStsTrack*) fStsTrkArray->UncheckedAt(sts_trk_idx))->GetNofMvdHits()
+ : 0;
+
+ int32_t sts_trk_size = sts_trk_idx != -1 && fStsTrkArray != nullptr
+ ? ((CbmStsTrack*) fStsTrkArray->UncheckedAt(sts_trk_idx))->GetNofStsHits()
+ : 0;
+
+ int32_t rich_trk_size = rich_trk_idx != -1 && fRchTrkArray != nullptr
+ ? ((CbmTrack*) fRchTrkArray->UncheckedAt(rich_trk_idx))->GetNofHits()
+ : 0;
+
+ int32_t much_trk_size = much_trk_idx != -1 && fMchTrkArray != nullptr
+ ? ((CbmTrack*) fMchTrkArray->UncheckedAt(much_trk_idx))->GetNofHits()
+ : 0;
+
+ int32_t trd_trk_size = trd_trk_idx != -1 && fTrdTrkArray != nullptr
+ ? ((CbmTrack*) fTrdTrkArray->UncheckedAt(trd_trk_idx))->GetNofHits()
+ : 0;
+
+ int32_t tof_trk_size = tof_trk_idx != -1 && fTofTrkArray != nullptr
+ ? ((CbmTrack*) fTofTrkArray->UncheckedAt(tof_trk_idx))->GetNofHits()
+ : 0;
+
+ if (fAnalysisCuts != nullptr
+ && (!fAnalysisCuts->Check(CbmCutId::kGlobalTrackMvdSize, mvd_trk_size)
+ || !fAnalysisCuts->Check(CbmCutId::kGlobalTrackStsSize, sts_trk_size)
+ || !fAnalysisCuts->Check(CbmCutId::kGlobalTrackRichSize, rich_trk_size)
+ || !fAnalysisCuts->Check(CbmCutId::kGlobalTrackMuchSize, much_trk_size)
+ || !fAnalysisCuts->Check(CbmCutId::kGlobalTrackTrdSize, trd_trk_size)
+ || !fAnalysisCuts->Check(CbmCutId::kGlobalTrackTofSize, tof_trk_size)
+ || !fAnalysisCuts->Check(CbmCutId::kGlobalTrackChi2, trk->GetChi2())
+ || !fAnalysisCuts->Check(CbmCutId::kGlobalTrackPval, trk_p_val))) {
+ return;
+ }
+ LOG(debug) << Form("Processing %d StsHit that were attached to a CATrack", sts_trk_size);
+ CbmStsTrack* sts_track = (CbmStsTrack*) fStsTrkArray->UncheckedAt(sts_trk_idx);
+
+ // ------------------------------------
+ // Process StsHit attached to a CATrack
+ // ------------------------------------
+ for (int hit_idx = 0; hit_idx < sts_trk_size; hit_idx++) {
+ uint32_t sts_hit_idx = sts_track->GetStsHitIndex(hit_idx);
+ CbmStsHit* sts_hit = (CbmStsHit*) fStsHitArray->UncheckedAt(sts_hit_idx);
+ ProcessHit(sts_hit, true);
+ }
+ // ------------------------------------
+}
+
+void CbmStsHitAna::Exec(Option_t*)
+{
+ if (fCbmEvtArray != nullptr) {
+ uint32_t nb_events = fCbmEvtArray->GetEntriesFast();
+ double avg_nb_sts_hits = 0;
+ double avg_nb_of_glob_trk = 0;
+ for (uint32_t evt_idx = 0; evt_idx < nb_events; evt_idx++) {
+ ProcessEvent((CbmEvent*) fCbmEvtArray->UncheckedAt(evt_idx));
+ avg_nb_sts_hits += ((CbmEvent*) fCbmEvtArray->UncheckedAt(evt_idx))->GetNofData(ECbmDataType::kStsHit);
+ avg_nb_of_glob_trk += ((CbmEvent*) fCbmEvtArray->UncheckedAt(evt_idx))->GetNofData(ECbmDataType::kGlobalTrack);
+ }
+ LOG_IF(info, nb_events) << "Running CbmStsHitAna - Event-like - Entry: " << entry_ << "\tEvents: " << nb_events
+ << "\tAvg StsHit/Events: " << avg_nb_sts_hits / nb_events
+ << "\tAvg GlobalTrk/Events: " << avg_nb_of_glob_trk / nb_events;
+ LOG_IF(info, !nb_events) << "Running CbmStsHitAna - Event-like - Entry: " << entry_ << "\tEvents: " << nb_events;
+ }
+ else {
+ uint32_t n_of_hits = fStsHitArray->GetEntriesFast();
+ LOG(info) << "Running CbmStsHitAna - Time-like - Entry: " << entry_ << "\tStsHits: " << n_of_hits;
+ for (uint32_t hit_idx = 0; hit_idx < n_of_hits; hit_idx++) {
+ CbmStsHit* sts_hit = (CbmStsHit*) fStsHitArray->UncheckedAt(hit_idx);
+ ProcessHit(sts_hit, false);
+ }
+ }
+ entry_++;
+}
+
+void CbmStsHitAna::Finish()
+{
+ SaveToFile();
+ fH1D.clear();
+ fH2D.clear();
+}
diff --git a/analysis/detectors/sts/CbmStsHitAna.h b/analysis/detectors/sts/CbmStsHitAna.h
new file mode 100644
index 0000000000000000000000000000000000000000..fcefecd4e4c761a94325302f90e03e794411f0c6
--- /dev/null
+++ b/analysis/detectors/sts/CbmStsHitAna.h
@@ -0,0 +1,106 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+#ifndef CBMSTSHITANA_H
+#define CBMSTSHITANA_H
+
+#include "CbmEvent.h"
+#include "CbmGlobalTrack.h"
+#include "CbmStsAnaBase.h"
+#include "CbmStsCluster.h"
+#include "CbmStsHit.h"
+#include "CbmStsParSetModule.h"
+#include "CbmStsUtils.h"
+
+#include <FairTask.h>
+
+#include <TClonesArray.h>
+
+#include <cstring>
+#include <map>
+#include <unordered_set>
+
+/**
+ * @class CbmStsHitAna
+ * @brief Task for detailed charge-size analysis of STS hits.
+ *
+ * This class inherits from FairTask and CbmStsAnaBase and is designed to
+ * perform detailed analysis of charge and size correlations in STS (Silicon Tracking System) hits.
+ *
+ * Key functionalities include:
+ *
+ * - **Charge Correlation**:
+ * Analyzes the correlation of charge between the p-side and n-side of STS sensors.
+ *
+ * - **Time Correlation**:
+ * Analyzes the correlation of time between the p-side and n-side cluster forming a hit.
+ *
+ * - **Cluster Size Correlation**:
+ * Evaluates the size correlation between p-side and n-side clusters.
+ *
+ * - **Charge vs. Size Analysis**:
+ * Performs charge versus cluster size studies separately for each sensor side.
+ *
+ * - **Hit Type Separation**:
+ * All analyses are performed separately for:
+ * - **All hits** (all) - includes all reconstructed hits.
+ * - **Track hits** (trk) - includes only hits associated with reconstructed tracks.
+ */
+
+class CbmStsHitAna : public FairTask, public CbmStsAnaBase {
+ public:
+ CbmStsHitAna();
+ ~CbmStsHitAna() = default;
+
+ /** \brief Constructor
+ * \param cal_par_file charge calibration file
+ */
+ CbmStsHitAna(std::string cal_par_file);
+
+ InitStatus Init();
+ void Exec(Option_t*);
+ void Finish();
+
+ private:
+ uint32_t fMaxCluSize{10};
+ std::vector<const char*> fHitModifier;
+
+ std::unique_ptr<CbmStsParSetModule> fModuleParSet;
+ std::string fCalibrationFile;
+
+ TClonesArray* fCbmEvtArray{nullptr};
+
+ TClonesArray* fGlbTrkArray{nullptr};
+ TClonesArray* fStsTrkArray{nullptr};
+ TClonesArray* fRchTrkArray{nullptr};
+ TClonesArray* fMchTrkArray{nullptr};
+ TClonesArray* fTrdTrkArray{nullptr};
+ TClonesArray* fTofTrkArray{nullptr};
+
+ TClonesArray* fStsHitArray{nullptr};
+ TClonesArray* fStsCluArray{nullptr};
+
+ void BookHistograms(int32_t);
+
+ /** \brief Process an Cbm events
+ * It filters event based on the provided CbmCutMap
+ */
+ void ProcessEvent(CbmEvent*);
+
+ /** \brief Process an STS hit
+ * It filters hits based on the provided CbmCutMap and fill the corresponding histogram
+ * \param hit StsHit to proccess
+ * \param belong_to_trk Specify that the hit belongs to a track, branching the histogram filling
+ */
+ void ProcessHit(CbmStsHit* hit, bool belong_to_trk = false);
+
+ /** \brief Process an STS hit attached to a Track
+ * It filters hits based on the provided CbmCutMap
+ * and check the track filters
+ */
+ void ProcessGlobalTrack(CbmGlobalTrack*);
+
+ ClassDef(CbmStsHitAna, 1);
+};
+#endif
diff --git a/analysis/detectors/sts/CbmStsRecoBeamSpot.cxx b/analysis/detectors/sts/CbmStsRecoBeamSpot.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..16ed149009b6a87a30ed66b1a7d69cacca87561f
--- /dev/null
+++ b/analysis/detectors/sts/CbmStsRecoBeamSpot.cxx
@@ -0,0 +1,205 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+#include "CbmStsRecoBeamSpot.h"
+
+#include "CbmStsAddress.h"
+
+void CbmStsRecoBeamSpot::AddTarget(CbmTarget* trg) { AddTarget("", trg); }
+void CbmStsRecoBeamSpot::AddTarget(std::string trg_name, CbmTarget* trg)
+{
+ if (trg_name.length() == 0) {
+ trg_name = Form("target_%ld", fTargets.size());
+ }
+ fTargets[trg_name] = trg;
+}
+
+void CbmStsRecoBeamSpot::BookHistograms()
+{
+ std::string h_name;
+ std::string h_name_base = "";
+ for (auto& [element_i, geo_i] : fStsGeoInfo) {
+ if (element_i < 8) continue;
+
+ for (auto& [element_j, geo_j] : fStsGeoInfo) {
+ if (element_j < 8) continue;
+ int32_t unit_i = CbmStsAddress::GetElementId(element_i, kStsUnit);
+ int32_t unit_j = CbmStsAddress::GetElementId(element_j, kStsUnit);
+
+ if (unit_i >= unit_j) continue;
+
+ for (auto& [trg_name, trg] : fTargets) {
+ h_name_base = Form("0x%x:0x%x_%s", element_i, element_j, trg_name.c_str());
+
+ fSampleRangeXmin = trg->GetPosition().Px() - 10.; // HARDCODE
+ fSampleRangeXmax = trg->GetPosition().Px() + 10.; // HARDCODE
+ fSampleRangeYmin = trg->GetPosition().Py() - 10.; // HARDCODE
+ fSampleRangeYmax = trg->GetPosition().Py() + 10.; // HARDCODE
+ fSampleRangeZmin = trg->GetPosition().Pz() - 10.; // HARDCODE
+ fSampleRangeZmax = trg->GetPosition().Pz() + 10.; // HARDCODE
+
+
+ unsigned int nb_bins_x = (fSampleRangeXmax - fSampleRangeXmin) / fSampleBinSizeX;
+ unsigned int nb_bins_y = (fSampleRangeYmax - fSampleRangeYmin) / fSampleBinSizeY;
+ unsigned int nb_bins_z = (fSampleRangeZmax - fSampleRangeZmin) / fSampleBinSizeZ;
+
+ h_name = Form("%s_Y_vs_X", h_name_base.c_str());
+ fH2D[h_name] = std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), nb_bins_x, fSampleRangeXmin,
+ fSampleRangeXmax, nb_bins_y, fSampleRangeYmin, fSampleRangeYmax);
+ fH2D[h_name]->GetXaxis()->SetTitle("X [cm]");
+ fH2D[h_name]->GetYaxis()->SetTitle("Y [cm]");
+
+ h_name = Form("%s_X_vs_Z", h_name_base.c_str());
+ fH2D[h_name] = std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), nb_bins_z, fSampleRangeZmin,
+ fSampleRangeZmax, nb_bins_x, fSampleRangeXmin, fSampleRangeXmax);
+ fH2D[h_name]->GetXaxis()->SetTitle("Z [cm]");
+ fH2D[h_name]->GetYaxis()->SetTitle("X [cm]");
+
+ h_name = Form("%s_Y_vs_Z", h_name_base.c_str());
+ fH2D[h_name] = std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), nb_bins_z, fSampleRangeZmin,
+ fSampleRangeZmax, nb_bins_y, fSampleRangeYmin, fSampleRangeYmax);
+ fH2D[h_name]->GetXaxis()->SetTitle("Z [cm]");
+ fH2D[h_name]->GetYaxis()->SetTitle("Y [cm]");
+
+ h_name = Form("%s_Y_beam_vs_X_beam", h_name_base.c_str());
+ fH2D[h_name] = std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), nb_bins_x, fSampleRangeXmin,
+ fSampleRangeXmax, nb_bins_y, fSampleRangeYmin, fSampleRangeYmax);
+ fH2D[h_name]->GetXaxis()->SetTitle("X_beam [cm]");
+ fH2D[h_name]->GetYaxis()->SetTitle("Y_beam [cm]");
+
+ } // end targets loop
+ } // end sts_i modules loop
+ } // end sts_j modules loop
+}
+
+TVector3 CbmStsRecoBeamSpot::ExtrapolateTrackTo(CbmPixelHit* hit_0, CbmPixelHit* hit_1, CbmTarget* trg)
+{
+ TVector3 n0 = trg->GetPosition();
+ TVector3 n1(.0, .0, 1.);
+ n1.RotateY(trg->GetRotation());
+
+ TVector3 l0(hit_0->GetX(), hit_0->GetY(), hit_0->GetZ());
+ TVector3 l1(hit_1->GetX() - hit_0->GetX(), hit_1->GetY() - hit_0->GetY(), hit_1->GetZ() - hit_0->GetZ());
+
+ double l1_dot_n1 = l1.Dot(n1);
+ if (l1_dot_n1 <= 1e-6) { // tracklet is parallel (or close to it) to target plane
+ throw std::invalid_argument("Track-let parallel to target plane");
+ }
+
+ double t = (n0 - l0).Dot(n1) / l1.Dot(n1);
+ return l0 + t * l1;
+}
+
+void CbmStsRecoBeamSpot::BeamSpotReco()
+{
+ if (!fStsHits.size()) return;
+ for (int sts_unit_i = 0; sts_unit_i < nb_sts_station_ - 1; sts_unit_i++) {
+ for (int sts_unit_j = sts_unit_i + 1; sts_unit_j < nb_sts_station_; sts_unit_j++) {
+ for (auto sts_hit_i : fStsHits[sts_unit_i]) {
+ for (auto sts_hit_j : fStsHits[sts_unit_j]) {
+ int32_t sts_address_i = sts_hit_i->GetAddress();
+ int32_t sts_address_j = sts_hit_j->GetAddress();
+
+ for (auto& [trg_name, trg] : fTargets) {
+ TVector3 sol = ExtrapolateTrackTo(sts_hit_i, sts_hit_j, trg);
+ std::string h_name_base = Form("0x%x:0x%x_%s", sts_address_i, sts_address_j, trg_name.c_str());
+
+ std::string h_name = Form("%s_X_vs_Z", h_name_base.c_str());
+ fH2D[h_name]->Fill(sol.Pz(), sol.Px());
+
+ h_name = Form("%s_Y_vs_Z", h_name_base.c_str());
+ fH2D[h_name]->Fill(sol.Pz(), sol.Py());
+
+ h_name = Form("%s_Y_vs_X", h_name_base.c_str());
+ fH2D[h_name]->Fill(sol.Px(), sol.Py());
+
+ TVector3 sol_beam = sol - trg->GetPosition();
+
+ double x_beam = sol_beam.Px() / cos(trg->GetRotation());
+ double y_beam = sol_beam.Py();
+ h_name = Form("%s_Y_beam_vs_X_beam", h_name_base.c_str());
+ fH2D[h_name]->Fill(x_beam, y_beam);
+ }
+ }
+ }
+ }
+ }
+ fStsHits.clear();
+}
+
+
+void CbmStsRecoBeamSpot::FinishTask() { SaveToFile(); }
+
+
+void CbmStsRecoBeamSpot::ProcessEvent(CbmEvent* evt)
+{
+ int nb_sts_hits = evt->GetNofData(ECbmDataType::kStsHit);
+ if (nb_sts_hits > 0)
+ for (int hit_evt_idx = 0; hit_evt_idx < nb_sts_hits; hit_evt_idx++) {
+ int sts_hit_idx = evt->GetIndex(ECbmDataType::kStsHit, hit_evt_idx);
+ CbmStsHit* hit = (CbmStsHit*) fStsHitArray->UncheckedAt(sts_hit_idx);
+ ProcessHit(hit);
+ }
+}
+
+void CbmStsRecoBeamSpot::ProcessHit(CbmStsHit* hit)
+{
+ if (hit == nullptr) return;
+ int32_t address = hit->GetAddress();
+ int32_t unit = CbmStsAddress::GetElementId(address, kStsUnit);
+
+ if (fAnalysisCuts != nullptr && fAnalysisCuts->CheckStsHit(hit, fStsCluArray)) {
+ fStsHits[unit].push_back(hit);
+ }
+}
+
+void CbmStsRecoBeamSpot::Exec(Option_t*)
+{
+
+ // Check if CbmEvent
+ if (fCbmEvtArray != nullptr) {
+
+ uint32_t nb_events = fCbmEvtArray->GetEntriesFast();
+ LOG(info) << "Running CbmStsRecoBeamSpot - Event like - Entry:" << entry_ << "\tEvents: " << nb_events;
+ for (uint32_t evt_idx = 0; evt_idx < nb_events; evt_idx++) {
+ CbmEvent* event = (CbmEvent*) fCbmEvtArray->UncheckedAt(evt_idx);
+ if (fAnalysisCuts != nullptr && fAnalysisCuts->CheckEvent(event)) {
+ ProcessEvent(event);
+ BeamSpotReco();
+ }
+ }
+ }
+ else {
+
+ uint32_t n_of_hits = fStsHitArray->GetEntriesFast();
+ LOG(info) << "Running CbmStsRecoBeamSpot - Time like - Entry:" << entry_ << "\tStsHit: " << n_of_hits;
+ for (uint32_t hit_idx = 0; hit_idx < n_of_hits; hit_idx++) {
+ CbmStsHit* sts_hit = (CbmStsHit*) fStsHitArray->UncheckedAt(hit_idx);
+ ProcessHit(sts_hit);
+ }
+ BeamSpotReco();
+ }
+ entry_++;
+}
+
+InitStatus CbmStsRecoBeamSpot::Init()
+{
+ LOG(debug) << "Init CbmStsRecoBeamSpot ...";
+
+ FairRootManager* ioman = FairRootManager::Instance();
+ if (ioman == nullptr) return kERROR;
+
+ fCbmEvtArray = (TClonesArray*) ioman->GetObject("CbmEvent");
+ fStsHitArray = (TClonesArray*) ioman->GetObject("StsHit");
+ fStsCluArray = (TClonesArray*) ioman->GetObject("StsCluster");
+
+ for (auto& [trg_name, trg] : fTargets) {
+ LOG(debug) << trg_name << ": " << trg->ToString();
+ }
+
+ LoadSetup();
+ BookHistograms();
+
+ return kSUCCESS;
+}
diff --git a/analysis/detectors/sts/CbmStsRecoBeamSpot.h b/analysis/detectors/sts/CbmStsRecoBeamSpot.h
new file mode 100644
index 0000000000000000000000000000000000000000..8ea1172e795931794beb0e8af197e16d01a08ab0
--- /dev/null
+++ b/analysis/detectors/sts/CbmStsRecoBeamSpot.h
@@ -0,0 +1,113 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+#ifndef CBMSTSRECOBEAMSPOT_H
+#define CBMSTSRECOBEAMSPOT_H
+
+#include "CbmEvent.h"
+#include "CbmGlobalTrack.h"
+#include "CbmStsAnaBase.h"
+#include "CbmStsCluster.h"
+#include "CbmStsHit.h"
+#include "CbmStsTrack.h"
+#include "CbmStsUtils.h"
+#include "CbmTarget.h"
+#include "CbmTrack.h"
+
+#include <FairTask.h>
+
+#include <TClonesArray.h>
+#include <TVector3.h>
+
+#include <cstring>
+#include <map>
+#include <unordered_map>
+#include <unordered_set>
+
+/** \class CbmStsRecoBeamSpot
+ * \brief Task for reconstructing the beam spot using STS hits.
+ *
+ * This class inherits from FairTask and CbmStsAnaBase.
+ * It provides functionality for reconstructing the beam spot at
+ * different planes defined by CbmTarget
+ * using STS hits.
+ */
+class CbmStsRecoBeamSpot : public FairTask, public CbmStsAnaBase {
+ public:
+ CbmStsRecoBeamSpot() = default;
+ ~CbmStsRecoBeamSpot() = default;
+
+ /**
+ * @brief Add a CbmTarget object to the list of targets.
+ *
+ * Position and rotation are the only relevant members of such class for this task.
+ *
+ * @param target CbmTarget ptr to be added
+ */
+
+ void AddTarget(CbmTarget* target = nullptr);
+
+ /**
+ * @brief Add a CbmTarget object to the list of targets with a key as trg_name.
+ *
+ * @param trg_name key under which the @param target will be stored.
+ * If the key is an empty string, added targets will be named as target_<n_of_targets>.
+ */
+ void AddTarget(std::string trg_name = "", CbmTarget* target = nullptr);
+
+ void Exec(Option_t*);
+ InitStatus Init();
+ void FinishTask();
+
+ private:
+ /* Beam Spot sampling range
+ * It is dynamically caculate when multiple targets planes are provided
+ */
+ double fSampleRangeXmin{-10}; // cm
+ double fSampleRangeXmax{+10}; // cm
+ double fSampleRangeYmin{-10}; // cm
+ double fSampleRangeYmax{+10}; // cm
+ double fSampleRangeZmin{-10}; // cm
+ double fSampleRangeZmax{+10}; // cm
+
+ /* Bin width for Beam Spot 2D histograms */
+ double fSampleBinSizeX{0.01};
+ double fSampleBinSizeY{0.01};
+ double fSampleBinSizeZ{0.01};
+
+ std::vector<CbmStsTrack*> fStsTrks;
+ std::map<int32_t, std::vector<CbmStsHit*>> fStsHits;
+
+ std::map<std::string, CbmTarget*> fTargets;
+
+
+ TClonesArray* fCbmEvtArray{nullptr};
+ TClonesArray* fStsTrkArray{nullptr};
+ TClonesArray* fStsHitArray{nullptr};
+ TClonesArray* fStsCluArray{nullptr};
+
+ /** \brief Reconstruct the beam spot at each target planes */
+ void BeamSpotReco();
+
+ void BookHistograms();
+
+ /** \brief Extrapolate a track-let to a target plane*/
+ TVector3 ExtrapolateTrackTo(CbmPixelHit*, CbmPixelHit*, CbmTarget*);
+
+ /** \brief Process an Cbm events
+ * It filters event based on the provided CbmCutMap
+ */
+ void ProcessEvent(CbmEvent*);
+
+ /** \brief Process an STS track */
+ void ProcessStsTrack(CbmGlobalTrack*);
+
+ /** \brief Process an STS hit
+ * It filters hits based on the provided CbmCutMap
+ */
+ void ProcessHit(CbmStsHit*);
+
+ ClassDef(CbmStsRecoBeamSpot, 1);
+};
+#endif
diff --git a/analysis/detectors/sts/CbmStsResolution.cxx b/analysis/detectors/sts/CbmStsResolution.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..ba67102876af191051be05f5bf1702e3a4ed0f1c
--- /dev/null
+++ b/analysis/detectors/sts/CbmStsResolution.cxx
@@ -0,0 +1,355 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+#include "CbmStsResolution.h"
+
+#include "CbmPixelHit.h"
+#include "CbmStsAddress.h"
+#include "CbmStsTrack.h"
+#include "TPaveStats.h"
+#include "TStyle.h"
+
+CbmStsResolution::CbmStsResolution(double d_max) : fImpactParMax(d_max)
+{
+ LOG(debug) << "Creating an instance of CbmStsResolution ...";
+}
+
+void CbmStsResolution::BookHistograms()
+{
+ for (auto& [element, geo] : fStsGeoInfo) {
+ std::string h_name_base = element < 8 ? Form("Sts%d", element) : Form("Sts0x%x", element);
+
+ LOG(debug) << Form("Booking for %s", h_name_base.c_str());
+
+ double x_min = geo[0];
+ double x_max = geo[1];
+ double y_min = geo[2];
+ double y_max = geo[3];
+
+ unsigned int nb_bins_x = (x_max - x_min) / cbm_sts_utils::kStsDx;
+ unsigned int nb_bins_y = (y_max - y_min) / cbm_sts_utils::kStsDy;
+
+
+ double dr = 5e-4;
+ double r_min = -2.7 - 0.5 * dr;
+ double r_max = +2.7 + 0.5 * dr;
+ unsigned int nb_bins_r = (r_max - r_min) / dr;
+
+ for (auto mod : {":all", ":trk"}) {
+ std::string h_name;
+ h_name = Form("%s_dx_vs_x%s", h_name_base.c_str(), mod);
+ fH2D[h_name] = std::make_unique<TH2D>(h_name.c_str(), "", nb_bins_r, r_min, r_max, nb_bins_x, x_min, x_max);
+ fH2D[h_name]->GetXaxis()->SetTitle(Form("X_{trk} - X_{hit} [cm]"));
+ fH2D[h_name]->GetYaxis()->SetTitle("X_{hit} [cm]");
+ fH2D[h_name]->GetXaxis()->CenterTitle(true);
+ fH2D[h_name]->GetYaxis()->CenterTitle(true);
+
+ h_name = Form("%s_dx_vs_y%s", h_name_base.c_str(), mod);
+ fH2D[h_name] = std::make_unique<TH2D>(h_name.c_str(), "", nb_bins_r, r_min, r_max, nb_bins_y, y_min, y_max);
+ fH2D[h_name]->GetXaxis()->SetTitle(Form("X_{trk} - X_{hit} [cm]"));
+ fH2D[h_name]->GetYaxis()->SetTitle("Y_{hit} [cm]");
+ fH2D[h_name]->GetXaxis()->CenterTitle(true);
+ fH2D[h_name]->GetYaxis()->CenterTitle(true);
+
+ h_name = Form("%s_dy_vs_x%s", h_name_base.c_str(), mod);
+ fH2D[h_name] = std::make_unique<TH2D>(h_name.c_str(), "", nb_bins_r, r_min, r_max, nb_bins_x, x_min, x_max);
+ fH2D[h_name]->GetXaxis()->SetTitle(Form("Y_{trk} - Y_{hit} [cm]"));
+ fH2D[h_name]->GetYaxis()->SetTitle("X_{hit} [cm]");
+ fH2D[h_name]->GetXaxis()->CenterTitle(true);
+ fH2D[h_name]->GetYaxis()->CenterTitle(true);
+
+ h_name = Form("%s_dy_vs_y%s", h_name_base.c_str(), mod);
+ fH2D[h_name] = std::make_unique<TH2D>(h_name.c_str(), "", nb_bins_r, r_min, r_max, nb_bins_y, y_min, y_max);
+ fH2D[h_name]->GetXaxis()->SetTitle(Form("Y_{trk} - Y_{hit} [cm]"));
+ fH2D[h_name]->GetYaxis()->SetTitle("Y_{hit} [cm]");
+ fH2D[h_name]->GetXaxis()->CenterTitle(true);
+ fH2D[h_name]->GetYaxis()->CenterTitle(true);
+ }
+ if (element < 8) {
+ std::string h_name = Form("Sts%d_ref_hits", element);
+ fH2D[h_name] = std::make_unique<TH2D>(h_name.c_str(), "", nb_bins_x, x_min, x_max, nb_bins_y, y_min, y_max);
+ fH2D[h_name]->GetXaxis()->SetTitle(Form("X_{TrkHit} [cm]"));
+ fH2D[h_name]->GetYaxis()->SetTitle("Y_{TrkHit} [cm]");
+ fH2D[h_name]->GetXaxis()->CenterTitle(true);
+ fH2D[h_name]->GetYaxis()->CenterTitle(true);
+ }
+ }
+}
+
+void CbmStsResolution::ProcessEvent(CbmEvent* evt)
+{
+ if (fAnalysisCuts != nullptr && !fAnalysisCuts->CheckEvent(evt)) return;
+
+ int nb_sts_hits = evt->GetNofData(ECbmDataType::kStsHit);
+ LOG(debug) << Form("Processing event with %d StsHit", nb_sts_hits);
+ for (int hit_evt_idx = 0; hit_evt_idx < nb_sts_hits; hit_evt_idx++) {
+ int sts_hit_idx = evt->GetIndex(ECbmDataType::kStsHit, hit_evt_idx);
+ CbmStsHit* hit = (CbmStsHit*) fStsHitArray->UncheckedAt(sts_hit_idx);
+ ProcessHit(hit);
+ }
+ int nb_of_glob_trk = evt->GetNofData(ECbmDataType::kGlobalTrack);
+ for (int evt_glob_trk_idx = 0; evt_glob_trk_idx < nb_of_glob_trk; evt_glob_trk_idx++) {
+ int glob_trk_idx = evt->GetIndex(ECbmDataType::kGlobalTrack, evt_glob_trk_idx);
+ CbmGlobalTrack* glob_trk = (CbmGlobalTrack*) fGlbTrkArray->UncheckedAt(glob_trk_idx);
+ ProcessGlobalTrack(glob_trk);
+ }
+ BuildResidual();
+}
+
+void CbmStsResolution::BuildResidual()
+{
+ int n_of_units = fStsHits.size();
+
+ for (auto trk : fGlbTrks) {
+ auto sts_trk_idx = trk->GetStsTrackIndex();
+ if (sts_trk_idx == -1) continue;
+
+ CbmStsTrack* sts_track = (CbmStsTrack*) fStsTrkArray->UncheckedAt(sts_trk_idx);
+ if (((CbmStsTrack*) fStsTrkArray->UncheckedAt(sts_trk_idx))->GetNofStsHits() != 2) continue;
+
+ CbmStsHit* hit_i = (CbmStsHit*) fStsHitArray->UncheckedAt(sts_track->GetStsHitIndex(0));
+ CbmStsHit* hit_j = (CbmStsHit*) fStsHitArray->UncheckedAt(sts_track->GetStsHitIndex(1));
+
+ // Charge cut to Sts track hits
+ if (fTrkHitQmin > 0
+ && (cbm_sts_utils::GetHitCharge(hit_i, fStsCluArray) < fTrkHitQmin
+ || cbm_sts_utils::GetHitCharge(hit_j, fStsCluArray) < fTrkHitQmin))
+ continue;
+
+ int unit_i = CbmStsAddress::GetElementId(hit_i->GetAddress(), kStsUnit);
+ int unit_j = CbmStsAddress::GetElementId(hit_j->GetAddress(), kStsUnit);
+
+ // Checkk track crossed eneabled sensors
+ if (fActiveRefSensors.size()) {
+ if (std::find(fActiveRefSensors.begin(), fActiveRefSensors.end(), hit_i->GetAddress()) == fActiveRefSensors.end())
+ continue;
+ if (std::find(fActiveRefSensors.begin(), fActiveRefSensors.end(), hit_j->GetAddress()) == fActiveRefSensors.end())
+ continue;
+ }
+
+ double x_i = hit_i->GetX();
+ double y_i = hit_i->GetY();
+ double z_i = hit_i->GetZ();
+
+ double x_j = hit_j->GetX();
+ double y_j = hit_j->GetY();
+ double z_j = hit_j->GetZ();
+
+ fH2D[Form("Sts%d_ref_hits", unit_i)]->Fill(x_i, y_i);
+ fH2D[Form("Sts%d_ref_hits", unit_j)]->Fill(x_j, y_j);
+
+ double x_0, y_0, t_x, t_y;
+ if (tracking_use_two_sts_) {
+ t_x = (x_i - x_j) / (z_i - z_j);
+ t_y = (y_i - y_j) / (z_i - z_j);
+ x_0 = x_i - t_x * z_i;
+ y_0 = y_i - t_y * z_i;
+ }
+ else {
+ t_x = trk->GetParamFirst()->GetTx();
+ t_y = trk->GetParamFirst()->GetTy();
+ x_0 = x_i - t_x * z_i;
+ y_0 = y_i - t_y * z_i;
+ }
+
+ double d_trk_trg = TMath::Sqrt(x_0 * x_0 + y_0 * y_0);
+
+ if (fImpactParMax > 0 && fImpactParMax < d_trk_trg) continue;
+
+ for (int uut = 0; uut < n_of_units; uut++) {
+ if (uut == unit_i || uut == unit_j) continue;
+
+ CbmStsHit* closest_hit = nullptr;
+ double min_dist = 1e+6;
+ for (auto hit_dut : fStsHits[uut]) {
+ double x_dut = hit_dut->GetX();
+ double y_dut = hit_dut->GetY();
+ double z_dut = hit_dut->GetZ();
+
+ double x_trk = x_i + t_x * (z_dut - z_i);
+ double y_trk = y_i + t_y * (z_dut - z_i);
+
+ double dx = x_trk - x_dut;
+ double dy = y_trk - y_dut;
+
+ fH2D[Form("Sts0x%x_dx_vs_x:all", hit_dut->GetAddress())]->Fill(dx, x_dut);
+ fH2D[Form("Sts0x%x_dx_vs_y:all", hit_dut->GetAddress())]->Fill(dx, y_dut);
+ fH2D[Form("Sts0x%x_dy_vs_x:all", hit_dut->GetAddress())]->Fill(dy, x_dut);
+ fH2D[Form("Sts0x%x_dy_vs_y:all", hit_dut->GetAddress())]->Fill(dy, y_dut);
+
+ fH2D[Form("Sts%d_dx_vs_x:all", uut)]->Fill(dx, x_dut);
+ fH2D[Form("Sts%d_dx_vs_y:all", uut)]->Fill(dx, y_dut);
+ fH2D[Form("Sts%d_dy_vs_x:all", uut)]->Fill(dy, x_dut);
+ fH2D[Form("Sts%d_dy_vs_y:all", uut)]->Fill(dy, y_dut);
+
+ if (min_dist > dx * dx + dy * dy) {
+ min_dist = dx * dx + dy * dy;
+ closest_hit = hit_dut;
+ }
+ } // hit_dut loop
+
+ if (closest_hit) {
+ double x_dut = closest_hit->GetX();
+ double y_dut = closest_hit->GetY();
+ double z_dut = closest_hit->GetZ();
+
+ double x_trk = x_i + t_x * (z_dut - z_i);
+ double y_trk = y_i + t_y * (z_dut - z_i);
+
+ double dx = x_trk - x_dut;
+ double dy = y_trk - y_dut;
+
+ fH2D[Form("Sts0x%x_dx_vs_x:trk", closest_hit->GetAddress())]->Fill(dx, x_dut);
+ fH2D[Form("Sts0x%x_dx_vs_y:trk", closest_hit->GetAddress())]->Fill(dx, y_dut);
+ fH2D[Form("Sts0x%x_dy_vs_x:trk", closest_hit->GetAddress())]->Fill(dy, x_dut);
+ fH2D[Form("Sts0x%x_dy_vs_y:trk", closest_hit->GetAddress())]->Fill(dy, y_dut);
+
+ fH2D[Form("Sts%d_dx_vs_x:trk", uut)]->Fill(dx, x_dut);
+ fH2D[Form("Sts%d_dx_vs_y:trk", uut)]->Fill(dx, y_dut);
+ fH2D[Form("Sts%d_dy_vs_x:trk", uut)]->Fill(dy, x_dut);
+ fH2D[Form("Sts%d_dy_vs_y:trk", uut)]->Fill(dy, y_dut);
+ }
+ } // uut loop
+ } // track loop
+
+ fStsHits.clear();
+ fGlbTrks.clear();
+}
+
+void CbmStsResolution::ProcessHit(CbmStsHit* hit)
+{
+ if (hit == nullptr) return;
+
+ if (fAnalysisCuts != nullptr && !fAnalysisCuts->CheckStsHit(hit, fStsCluArray)) return;
+
+ int32_t address = hit->GetAddress();
+ int32_t unit = CbmStsAddress::GetElementId(address, kStsUnit);
+
+ fStsHits[unit].push_back(hit);
+}
+
+void CbmStsResolution::EnableRefSensor(int32_t address) { fActiveRefSensors.push_back(address); }
+
+void CbmStsResolution::ProcessGlobalTrack(CbmGlobalTrack* trk)
+{
+ if (trk == nullptr) return;
+
+ auto sts_trk_idx = trk->GetStsTrackIndex();
+ auto rich_trk_idx = trk->GetRichRingIndex();
+ auto much_trk_idx = trk->GetMuchTrackIndex();
+ auto trd_trk_idx = trk->GetTrdTrackIndex();
+ auto tof_trk_idx = trk->GetTofTrackIndex();
+ float trk_p_val = TMath::Prob(trk->GetChi2(), trk->GetNDF());
+
+ // Apply GlobalTracks cuts
+ int32_t mvd_trk_size = sts_trk_idx != -1 && fStsTrkArray != nullptr
+ ? ((CbmStsTrack*) fStsTrkArray->UncheckedAt(sts_trk_idx))->GetNofMvdHits()
+ : 0;
+
+ int32_t sts_trk_size = sts_trk_idx != -1 && fStsTrkArray != nullptr
+ ? ((CbmStsTrack*) fStsTrkArray->UncheckedAt(sts_trk_idx))->GetNofStsHits()
+ : 0;
+
+ int32_t rich_trk_size = rich_trk_idx != -1 && fRchTrkArray != nullptr
+ ? ((CbmTrack*) fRchTrkArray->UncheckedAt(rich_trk_idx))->GetNofHits()
+ : 0;
+
+ int32_t much_trk_size = much_trk_idx != -1 && fMchTrkArray != nullptr
+ ? ((CbmTrack*) fMchTrkArray->UncheckedAt(much_trk_idx))->GetNofHits()
+ : 0;
+
+ int32_t trd_trk_size = trd_trk_idx != -1 && fTrdTrkArray != nullptr
+ ? ((CbmTrack*) fTrdTrkArray->UncheckedAt(trd_trk_idx))->GetNofHits()
+ : 0;
+
+ int32_t tof_trk_size = tof_trk_idx != -1 && fTofTrkArray != nullptr
+ ? ((CbmTrack*) fTofTrkArray->UncheckedAt(tof_trk_idx))->GetNofHits()
+ : 0;
+
+ if (fAnalysisCuts != nullptr
+ && (!fAnalysisCuts->Check(CbmCutId::kGlobalTrackMvdSize, mvd_trk_size)
+ || !fAnalysisCuts->Check(CbmCutId::kGlobalTrackStsSize, sts_trk_size)
+ || !fAnalysisCuts->Check(CbmCutId::kGlobalTrackRichSize, rich_trk_size)
+ || !fAnalysisCuts->Check(CbmCutId::kGlobalTrackMuchSize, much_trk_size)
+ || !fAnalysisCuts->Check(CbmCutId::kGlobalTrackTrdSize, trd_trk_size)
+ || !fAnalysisCuts->Check(CbmCutId::kGlobalTrackTofSize, tof_trk_size)
+ || !fAnalysisCuts->Check(CbmCutId::kGlobalTrackChi2, trk->GetChi2())
+ || !fAnalysisCuts->Check(CbmCutId::kGlobalTrackPval, trk_p_val))) {
+ return;
+ }
+
+ fGlbTrks.push_back(trk);
+}
+
+void CbmStsResolution::Exec(Option_t*)
+{
+ switch (fInputType) {
+ case kEventMode: {
+ if (fCbmEvtArray == nullptr) {
+ LOG(error) << "Invalid branch for event-mode: Branch CbmEvent -> nullptr";
+ break;
+ }
+
+ uint32_t nb_events = fCbmEvtArray->GetEntriesFast();
+ LOG(info) << "Processing entry " << entry_ << "\tCbmEvent: " << nb_events;
+
+ for (uint32_t evt_idx = 0; evt_idx < nb_events; evt_idx++) {
+ CbmEvent* event = (CbmEvent*) fCbmEvtArray->UncheckedAt(evt_idx);
+ ProcessEvent(event);
+ }
+ break;
+ }
+ case kMCTimeMode: {
+ break;
+ }
+ case kMCEventMode:
+ case kTimeMode:
+ default: {
+ uint32_t n_of_hits = fStsHitArray->GetEntriesFast();
+ uint32_t n_of_glb_trk = fGlbTrkArray->GetEntriesFast();
+ LOG(info) << "Processing entry " << entry_ << "\tStsHit: " << n_of_hits << "\tGlobalTrack: " << n_of_glb_trk;
+
+ for (uint32_t hit_idx = 0; hit_idx < n_of_hits; hit_idx++) {
+ ProcessHit((CbmStsHit*) fStsHitArray->UncheckedAt(hit_idx));
+ }
+
+ for (uint32_t glb_trk_idx = 0; glb_trk_idx < n_of_glb_trk; glb_trk_idx++) {
+ ProcessGlobalTrack((CbmGlobalTrack*) fGlbTrkArray->UncheckedAt(glb_trk_idx));
+ }
+
+ BuildResidual();
+ break;
+ }
+ }
+ entry_++;
+}
+
+InitStatus CbmStsResolution::Init()
+{
+ LOG(debug) << "Init CbmStsResolution ...";
+
+ FairRootManager* ioman = FairRootManager::Instance();
+ if (ioman != nullptr) {
+
+ fCbmEvtArray = (TClonesArray*) ioman->GetObject("CbmEvent");
+
+ fGlbTrkArray = (TClonesArray*) ioman->GetObject("GlobalTrack");
+ fStsTrkArray = (TClonesArray*) ioman->GetObject("StsTrack");
+ fRchTrkArray = (TClonesArray*) ioman->GetObject("RichTrack");
+ fMchTrkArray = (TClonesArray*) ioman->GetObject("MuchTrack");
+ fTrdTrkArray = (TClonesArray*) ioman->GetObject("TrdTrack");
+ fTofTrkArray = (TClonesArray*) ioman->GetObject("TofTrack");
+
+ fStsHitArray = (TClonesArray*) ioman->GetObject("StsHit");
+ fStsCluArray = (TClonesArray*) ioman->GetObject("StsCluster");
+ }
+
+ LoadSetup();
+
+ BookHistograms();
+
+ return kSUCCESS;
+}
+
+void CbmStsResolution::Finish() { SaveToFile(); }
diff --git a/analysis/detectors/sts/CbmStsResolution.h b/analysis/detectors/sts/CbmStsResolution.h
new file mode 100644
index 0000000000000000000000000000000000000000..a6938e855cc901c5db71e5617233f4cc7470427d
--- /dev/null
+++ b/analysis/detectors/sts/CbmStsResolution.h
@@ -0,0 +1,133 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+#ifndef CBMSTSRESOLUTION_H
+#define CBMSTSRESOLUTION_H
+
+#include "CbmEvent.h"
+#include "CbmGlobalTrack.h"
+#include "CbmStsAnaBase.h"
+#include "CbmStsCluster.h"
+#include "CbmStsHit.h"
+#include "CbmStsUtils.h"
+
+#include <FairTask.h>
+
+#include <TClonesArray.h>
+
+#include <cstring>
+#include <map>
+#include <unordered_map>
+#include <unordered_set>
+
+
+/**
+ * \class CbmStsResolution
+ * \brief Task for estimating the spatial resolution of the STS detector using track-hit residuals.
+ *
+ * The `CbmStsResolution` class is responsible for calculating the residuals between the hits
+ * in the Silicon Tracking System (STS) and the corresponding reconstructed tracks. These
+ * residuals are used to estimate the spatial resolution of the STS detector.
+ *
+ * ## Output:
+ * The class produces a set of 2D histograms showing the distribution of residuals
+ * as a function of the reconstructed hit coordinates. These histograms provide insight into
+ * the detector's resolution performance across different spatial regions.
+ *
+ * If functon EnableRefSensor() is called during analysis configuration, the residuals will be only produced
+ * when the track hit the selected sensors. This is especially useful to restrict bad performing sensors.
+ *
+ * ## Track Extrapolation Modes:
+ * The residuals can be computed using two approaches:
+ *
+ * 1. **Local STS-based extrapolation:** Uses only two STS hits to form a straight-line
+ * extrapolation to calculate the expected hit position.
+ *
+ * 2. **Global track-based extrapolation:** Utilizes the full set of global track parameters
+ * obtained from the reconstruction process to extrapolate the expected position at the hit.
+ *
+ * ## The residual are calculated in two modes:alignas
+ * - `":all"`: All hits are used for the residual calculation. (large combinatoric)
+ * - `":trk"`: Only the hits clsosest to the extrapolation point is used (selection bias, less combinatoric).
+ *
+ * ## Track selection:
+ * It can be configured using the `CbmCutMap` class.
+ * Additionally, the user can set a minimum charge threshold for the STS hits attached to the track
+ * and set a cut to the maximun value for the impact parameter.
+ *
+ * ## Future Extension:
+ * The class is planned to be extended with functionality to support **track refitting using
+ * only STS hits**, which would allow for improved and more localized resolution studies
+ * decoupled from the global tracking performance.
+ *
+ * Inherits from:
+ * - `FairTask`: Base class for FairRoot tasks, enabling integration with the FAIR framework.
+ * - `CbmStsAnaBase`: Provides STS-specific analysis utilities and interfaces.
+ *
+ */
+class CbmStsResolution : public FairTask, public CbmStsAnaBase {
+ public:
+ CbmStsResolution() = default;
+ ~CbmStsResolution() = default;
+ ;
+
+ /** \brief Constructor */
+ CbmStsResolution(double);
+
+ InitStatus Init();
+ void Exec(Option_t*);
+ void Finish();
+
+ void SetInputType(InputType type) { fInputType = type; };
+
+ void EnableRefSensor(int32_t);
+
+ void SetMinTrkHitQ(double q) { fTrkHitQmin = q; }
+
+ void TrackUseTwoSts() { tracking_use_two_sts_ = true; }
+
+ private:
+ double fImpactParMax{-1};
+ double fTrkHitQmin{-1};
+
+ std::vector<CbmGlobalTrack*> fGlbTrks;
+ std::map<int32_t, std::vector<CbmStsHit*>> fStsHits;
+
+ InputType fInputType = InputType::kEventMode;
+
+ bool tracking_use_two_sts_{false};
+
+ std::vector<int32_t> fActiveRefSensors;
+
+ TClonesArray* fCbmEvtArray{nullptr};
+
+ TClonesArray* fGlbTrkArray{nullptr};
+ TClonesArray* fStsTrkArray{nullptr};
+ TClonesArray* fRchTrkArray{nullptr};
+ TClonesArray* fMchTrkArray{nullptr};
+ TClonesArray* fTrdTrkArray{nullptr};
+ TClonesArray* fTofTrkArray{nullptr};
+
+ TClonesArray* fStsHitArray{nullptr};
+ TClonesArray* fStsCluArray{nullptr};
+
+ void BookHistograms();
+
+ /** \brief Process an Cbm events
+ * It filters event based on the provided CbmCutMap
+ */
+ void ProcessEvent(CbmEvent*);
+
+ /** \brief Process an STS hit
+ * It filters hits based on the provided CbmCutMap
+ */
+ void ProcessHit(CbmStsHit*);
+
+ void ProcessGlobalTrack(CbmGlobalTrack* trk);
+
+ void BuildResidual();
+
+ ClassDef(CbmStsResolution, 1);
+};
+#endif
diff --git a/analysis/detectors/sts/CbmStsTimeCal.cxx b/analysis/detectors/sts/CbmStsTimeCal.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..2df8ae30532f63219d0b9b976f2e387507bd4bbf
--- /dev/null
+++ b/analysis/detectors/sts/CbmStsTimeCal.cxx
@@ -0,0 +1,606 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+#include "CbmStsTimeCal.h"
+
+#include "CbmBmonDigi.h"
+#include "CbmFsdDigi.h"
+#include "CbmMuchDigi.h"
+#include "CbmMvdDigi.h"
+#include "CbmRichDigi.h"
+#include "CbmStsDigi.h"
+#include "CbmTofDigi.h"
+#include "CbmTrdDigi.h"
+
+#include <TColor.h>
+#include <TDirectory.h>
+#include <TF1.h>
+#include <TFitResult.h>
+#include <TLine.h>
+#include <TPaveText.h>
+#include <TStyle.h>
+
+#include <cassert>
+#include <ctime>
+#include <fstream>
+#include <iostream>
+#include <typeinfo>
+
+#include <fmt/core.h>
+#include <yaml-cpp/emitter.h>
+#include <yaml-cpp/emittermanip.h>
+#include <yaml-cpp/node/node.h>
+
+CbmStsTimeCal::CbmStsTimeCal(ECbmModuleId ref_sys, double t_min, double t_max)
+ : fTimeWindowMin(t_min)
+ , fTimeWindowMax(t_max)
+ , fRefSystem(ref_sys)
+{
+}
+
+CbmStsTimeCal::CbmStsTimeCal(int run_id, ECbmModuleId ref_sys, double t_min, double t_max)
+{
+ fRunId = run_id;
+ fTimeWindowMin = t_min;
+ fTimeWindowMax = t_max;
+ fRefSystem = ref_sys;
+}
+
+
+InitStatus CbmStsTimeCal::Init()
+{
+ FairRootManager* ioman = FairRootManager::Instance();
+ if (ioman != nullptr) {
+ fDigiManager = CbmDigiManager::Instance();
+ fDigiManager->Init();
+
+ fCbmEvtArray = (TClonesArray*) ioman->GetObject("CbmEvent");
+
+ if (!fDigiManager->IsPresent(ECbmModuleId::kSts)) LOG(fatal) << GetName() << ": No StsDigi branch in input!";
+ if (!fDigiManager->IsPresent(fRefSystem))
+ LOG(fatal) << GetName() << ": No " << ToString(fRefSystem) << " branch in input!";
+
+ LoadWalkFromFile();
+
+ fOutputPath = fRunId > 0 ? Form("%d", fRunId) : ".";
+ if (gSystem->AccessPathName(fOutputPath.c_str(), kFileExists)) {
+ if (system(("mkdir -p " + fOutputPath).c_str())) { // Check output folder
+ LOG(warning) << "Could not create output directory\n Setting output path at current location:\n";
+ }
+ }
+ fReportOffset = std::make_unique<TFile>(Form("%s/cbm_sts_time_cal_offset.root", fOutputPath.c_str()), "RECREATE");
+ fReportFile = std::make_unique<TFile>(Form("%s/cbm_sts_time_cal_report.root", fOutputPath.c_str()), "RECREATE");
+ fReportFit = std::make_unique<TFile>(Form("%s/cbm_sts_time_cal_fits.root", fOutputPath.c_str()), "RECREATE");
+
+ return kSUCCESS;
+ }
+ return kERROR;
+}
+
+void CbmStsTimeCal::InitTimeWalkMap(int32_t address)
+{
+ for (unsigned int asic_idx = 0; asic_idx < 16; asic_idx++) {
+ fWalkMapRaw[address][asic_idx] = std::vector<double>(32, 0);
+ auto loaded_par = fWalkMap.Get(address, asic_idx);
+ if (loaded_par.size() == 31) {
+ for (int adc = 1; adc <= 31; adc++) {
+ fWalkMapRaw[address][asic_idx][adc] = loaded_par[adc - 1];
+ }
+ }
+ }
+}
+
+void CbmStsTimeCal::SetWalkFile(std::string f_name)
+{
+ LOG(debug) << Form("Setting user define time calibration file: %s", f_name.c_str());
+ fParFile = f_name;
+}
+
+void CbmStsTimeCal::LoadWalkFromFile()
+{
+ if (!fParFile.length()) return;
+
+ if (TString(fParFile.c_str()).EndsWith(".yaml") || TString(fParFile.c_str()).EndsWith(".yml")) {
+ LOG(info) << Form("Loading time calibration from parameter file: %s", fParFile.c_str());
+ fWalkMap = cbm::algo::yaml::ReadFromFile<cbm::algo::sts::WalkMap>(fParFile);
+ fGlobalTimeOffset = fWalkMap.GetSystemTimeOffset();
+ LOG(info) << "Current system offset: " << fGlobalTimeOffset << " ns";
+ }
+}
+
+void CbmStsTimeCal::BookHistograms(int32_t address)
+{
+ LOG(debug) << Form("Booking histograms for module: 0x%x", address);
+
+ std::string h_name, h_title;
+ int nb_time_bins = (std::abs(fTimeWindowMax - fTimeWindowMin) + kStsClock) / kStsClock;
+
+ int unit = CbmStsAddress::GetElementId(address, kStsUnit);
+ int ladd = CbmStsAddress::GetElementId(address, kStsLadder);
+ int modu = CbmStsAddress::GetElementId(address, kStsModule);
+ std::string smart_name = Form("U%d_L%d_M%d", unit, ladd, modu);
+
+ for (int asic_idx = 0; asic_idx < 16; asic_idx++) {
+ h_name = Form("0x%x_%02d", address, asic_idx);
+ h_title = Form("%s_%02d", smart_name.c_str(), asic_idx);
+
+ fH2D[h_name] =
+ std::make_unique<TH2D>(h_name.c_str(), h_title.c_str(), nb_time_bins, fTimeWindowMin - 0.5 * kStsClock,
+ fTimeWindowMax + 0.5 * kStsClock, 31, 1, 32);
+
+ fH2D[h_name]->GetXaxis()->SetTitle(Form("t_{%s} - t_{StsDigi} [ns]", ToString(fRefSystem).c_str()));
+ fH2D[h_name]->GetYaxis()->SetTitle("Charge_{StsDigi} [ADC]");
+ }
+
+ h_name = Form("0x%x_dt_vs_channel", address);
+ h_title = smart_name;
+ fH2D[h_name] = std::make_unique<TH2D>(h_name.c_str(), h_title.c_str(), 2048, 0, 2048, nb_time_bins,
+ fTimeWindowMin - 0.5 * kStsClock, fTimeWindowMax + 0.5 * kStsClock);
+
+ fH2D[h_name]->GetYaxis()->SetTitle(Form("t_{%s} - t_{StsDigi} [ns]", ToString(fRefSystem).c_str()));
+ fH2D[h_name]->GetXaxis()->SetTitle("Channel_{StsDigi}");
+
+ fAddressBook.insert(address);
+}
+
+void CbmStsTimeCal::Exec(Option_t*)
+{
+ LOG(info) << "Running CbmStsTimeCal - Entry " << entry_;
+
+ if (fRefSystem == ECbmModuleId::kBmon) CheckTiming<CbmBmonDigi>();
+ if (fRefSystem == ECbmModuleId::kTof) CheckTiming<CbmTofDigi>();
+ if (fRefSystem == ECbmModuleId::kRich) CheckTiming<CbmRichDigi>();
+
+ entry_++;
+}
+
+template<class Digi>
+void CbmStsTimeCal::CheckTiming()
+{
+ auto sts_digis_ = fDigiManager->GetArray<CbmStsDigi>();
+ auto ref_digis_ = fDigiManager->GetArray<Digi>();
+
+ size_t nb_sts_digis = fDigiManager->GetNofDigis(ECbmModuleId::kSts);
+ size_t nb_ref_digis = fDigiManager->GetNofDigis(Digi::GetSystem());
+
+ LOG(info) << nb_sts_digis << "\t" << nb_ref_digis;
+
+ // Loops over Ref Detector
+ for (size_t ref_digi_idx = 0; ref_digi_idx < nb_ref_digis - 1; ref_digi_idx++) {
+ const Digi* ref_digi = &ref_digis_[ref_digi_idx];
+ double ref_digi_time = ref_digi->GetTime();
+
+ if (fAnalysisCuts != nullptr
+ && !fAnalysisCuts->Check(CbmCutId::kBmonDigiSide, CbmTofAddress::GetChannelSide(ref_digi->GetAddress())))
+ continue;
+
+ double first_digi_in_window = ref_digi_time + fTimeWindowMin;
+ double last_digi_in_window = ref_digi_time + fTimeWindowMax;
+ // Find first_digi_in_window
+ int lo = 0, hi = nb_sts_digis - 1;
+ while (lo < hi) {
+ int m = (lo + hi) / 2;
+ if (sts_digis_[m].GetTime() >= first_digi_in_window) {
+ hi = m;
+ }
+ else {
+ lo = m + 1;
+ }
+ }
+ int lower_hitB = lo;
+
+ // Find last_digi_in_window
+ lo = 0, hi = nb_sts_digis - 1;
+ while (lo < hi) {
+ int m = (lo + hi + 1) / 2;
+ if (sts_digis_[m].GetTime() <= last_digi_in_window) {
+ lo = m;
+ }
+ else {
+ hi = m - 1;
+ }
+ }
+ int upper_hitB = lo;
+ if (lower_hitB > upper_hitB) {
+ LOG(debug) << "Not time match found for the current hit within the time window";
+ continue;
+ }
+
+ for (int sts_digi_idx = lower_hitB; sts_digi_idx < upper_hitB; sts_digi_idx++) {
+ const CbmStsDigi* sts_digi = &sts_digis_[sts_digi_idx];
+ int32_t sts_digi_addr = sts_digi->GetAddress();
+
+ if (!fAddressBook.count(sts_digi_addr)) {
+ BookHistograms(sts_digi_addr);
+ InitTimeWalkMap(sts_digi_addr);
+ }
+
+ int32_t sts_digi_chan = sts_digi->GetChannel();
+ int32_t sts_digi_char = sts_digi->GetCharge();
+ int32_t asic_idx = sts_digi_chan / 128;
+ double sts_digi_time = sts_digi->GetTime();
+
+ fH2D[Form("0x%x_%02d", sts_digi_addr, asic_idx)]->Fill(ref_digi_time - sts_digi_time, sts_digi_char);
+ fH2D[Form("0x%x_dt_vs_channel", sts_digi_addr)]->Fill(sts_digi_chan, ref_digi_time - sts_digi_time);
+ } // end sts loop
+ } // end ref loop
+}
+
+
+struct FitStsTime {
+ std::unique_ptr<TH1> h_x;
+ std::unique_ptr<TF1> f_x;
+ std::unique_ptr<TF1> f_peak;
+ std::unique_ptr<TF1> f_bkgd;
+ TFitResultPtr fit_ptr;
+ bool use_fit_result{false};
+ double time_offset{0};
+
+ TCanvas* Draw()
+ {
+ if (h_x == nullptr) {
+ return nullptr;
+ }
+ int pad_x = 1024;
+ int pad_y = 1024;
+
+ std::unique_ptr<TCanvas> canvas = std::make_unique<TCanvas>("_c0", "", pad_x, pad_y);
+ canvas->cd();
+
+ double line_width = 0.06;
+ double pave0_top_edge = 0.95;
+ double pave0_bot_edge = pave0_top_edge - line_width * 3;
+ double pave1_top_edge = pave0_bot_edge;
+ double pave1_bot_edge = pave1_top_edge - line_width * 3;
+
+ h_x->SetLineStyle(1);
+ h_x->SetLineWidth((0.002 * pad_y));
+ h_x->SetLineColor(kBlack);
+ h_x->DrawClone("histo");
+
+ if (use_fit_result) {
+ f_x->SetTitle("Fitted function");
+ f_x->SetLineStyle(1);
+ f_x->SetLineWidth((0.002 * pad_y));
+ f_x->SetLineColor(kRed);
+ f_x->DrawClone("same");
+ }
+
+ if (use_fit_result) {
+ f_peak->SetTitle("Peak");
+ f_peak->SetLineStyle(kDotted);
+ f_peak->SetLineWidth((0.002 * pad_y));
+ f_peak->SetLineColor(kBlue);
+ f_peak->DrawClone("same");
+ std::unique_ptr<TPaveText> pave = std::make_unique<TPaveText>(0.16, pave0_bot_edge, 0.45, pave0_top_edge, "NDC");
+ pave->SetBorderSize(1);
+ pave->SetFillColor(0);
+ pave->SetTextAlign(12);
+ pave->AddText(Form("p0 = %.3f", f_peak->GetParameter(0)))->SetTextColor(kBlue);
+ pave->AddText(Form("p1 = %.3f", f_peak->GetParameter(1)))->SetTextColor(kBlue);
+ pave->AddText(Form("p2 = %.3f", f_peak->GetParameter(2)))->SetTextColor(kBlue);
+ pave->DrawClone();
+ }
+
+ if (use_fit_result) {
+ f_bkgd->SetTitle("Background");
+ f_bkgd->SetLineStyle(kDotted);
+ f_bkgd->SetLineWidth((0.002 * pad_y));
+ f_bkgd->SetLineColor(kOrange);
+ f_bkgd->DrawClone("same");
+ std::unique_ptr<TPaveText> pave = std::make_unique<TPaveText>(0.16, pave1_bot_edge, 0.45, pave1_top_edge, "NDC");
+ pave->SetBorderSize(1);
+ pave->SetFillColor(0);
+ pave->SetTextAlign(12);
+ pave->AddText(Form("p0 = %.3f", f_bkgd->GetParameter(0)))->SetTextColor(kOrange);
+ pave->AddText(Form("p1 = %.3f", f_bkgd->GetParameter(1)))->SetTextColor(kOrange);
+ pave->AddText(Form("p2 = %.3f", f_bkgd->GetParameter(2)))->SetTextColor(kOrange);
+ pave->DrawClone();
+ }
+
+ TLine offset(time_offset, h_x->GetMinimum(), time_offset, h_x->GetMaximum());
+ offset.SetLineWidth(0.004 * pad_y);
+ offset.SetLineColor(kGray);
+ offset.SetLineStyle(kDotted);
+ offset.DrawClone();
+
+ // Fit convergence
+ std::unique_ptr<TPaveText> fit_pave = std::make_unique<TPaveText>(0.65, 0.75, 0.95, 0.95, "NDC");
+ fit_pave->SetBorderSize(1);
+ fit_pave->SetFillColor(0);
+ fit_pave->SetTextAlign(12);
+ fit_pave->AddText(Form("NDF = %d", fit_ptr->Ndf()));
+ fit_pave->AddText(Form("Chi2 = %.3f", fit_ptr->Chi2()));
+ fit_pave->AddText(Form("PRob = %E", fit_ptr->Prob()));
+ fit_pave->DrawClone();
+
+ return (TCanvas*) canvas->Clone();
+ }
+
+ FitStsTime(TH1D* h)
+ {
+ if (h == nullptr) {
+ return;
+ }
+ h_x = std::make_unique<TH1D>(*h);
+ h_x->Scale(1. / h_x->Integral());
+
+ // Fit configuration
+ double time_resolution = 4.8;
+ double h_max = h_x->GetBinCenter(h_x->GetMaximumBin());
+ double fit_x_min = h_max - 10. * time_resolution;
+ double fit_x_max = h_max + 10. * time_resolution;
+
+ f_x = std::make_unique<TF1>("f_x", "[0]*TMath::Gaus(x, [1], [2]) + [3] + [4]*x + [5]*x*x", fit_x_min, fit_x_max);
+ f_peak = std::make_unique<TF1>("f_peak", "[0]*TMath::Gaus(x, [1], [2])", fit_x_min, fit_x_max);
+ f_bkgd = std::make_unique<TF1>("f_bkgd", "[0] + [1]*x + [2]*x*x", fit_x_min, fit_x_max);
+
+ // Initial values
+ f_x->SetParameter(0, 0.5);
+ f_x->SetParameter(1, h_max);
+ f_x->SetParameter(2, time_resolution);
+ f_x->SetParameter(5, -1e-4);
+
+ // Limits
+ f_x->SetParLimits(0, 0.0, 1.0);
+ f_x->SetParLimits(1, h_max - 2. * time_resolution, h_max + 2. * time_resolution);
+ f_x->SetParLimits(2, 3.125, 25);
+ f_x->SetParLimits(5, -1, 0.0);
+
+
+ fit_ptr = h_x->Fit("f_x", "SQN0", "", fit_x_min, fit_x_max);
+
+ f_x->SetParameters(fit_ptr->GetParams());
+ f_peak->SetParameter(0, fit_ptr->Parameter(0));
+ f_peak->SetParameter(1, fit_ptr->Parameter(1));
+ f_peak->SetParameter(2, fit_ptr->Parameter(2));
+ f_bkgd->SetParameter(0, fit_ptr->Parameter(3));
+ f_bkgd->SetParameter(1, fit_ptr->Parameter(4));
+ f_bkgd->SetParameter(2, fit_ptr->Parameter(5));
+
+ bool par_out_of_limits =
+ fit_ptr->Parameter(1) < h_max - 2. * time_resolution || h_max + 2. * time_resolution < fit_ptr->Parameter(1);
+
+ use_fit_result = !fit_ptr && fit_ptr->IsValid() && !par_out_of_limits;
+
+ if (use_fit_result) {
+ time_offset = f_peak->GetParameter(1);
+ }
+ else {
+ time_offset = h_max;
+ }
+ }
+};
+
+double CbmStsTimeCal::FindModuleOffset(int32_t address)
+{
+
+ double module_offset = 0;
+ TH2D* dt_vs_chn = fH2D[Form("0x%x_dt_vs_channel", address)].get();
+ if (dt_vs_chn != nullptr) {
+ TH1D* p_y = (TH1D*) dt_vs_chn->ProjectionY();
+
+ double peak = p_y->GetBinCenter(p_y->GetMaximumBin());
+
+ FitStsTime offset_fit(p_y);
+ TCanvas* drawing = offset_fit.Draw();
+ if (drawing != nullptr) {
+ fReportOffset->cd();
+ drawing->Write(Form("0x%x_offset", address));
+ }
+
+ if (offset_fit.use_fit_result) {
+ module_offset = offset_fit.time_offset;
+ }
+ else {
+ module_offset = peak;
+ }
+ }
+ return module_offset;
+}
+
+void CbmStsTimeCal::CheckTimeWalk()
+{
+ /* This function can be optimized by parallelizing the loop over modules or ASIC idx, depending on the ratio
+ * For mCBM case, there are usually less than 16 modules so it makes more sense to parallelize the ASIC idx loop
+ * but for the full CBM case, the amount of modules greatly exceeds the number of ASIC per modules (16)
+ */
+ for (int32_t module : fAddressBook) { // loop over modules
+ double module_offset = FindModuleOffset(module);
+ for (unsigned int asic_idx = 0; asic_idx < 16; asic_idx++) { // loop over ASICs
+
+ double adc_channel[32];
+ double adc_channel_err[32];
+ double time_offset[32];
+ double time_offset_err[32];
+ int n_pts = 0;
+
+ std::string h_name = Form("0x%x_%02d", module, asic_idx);
+ if (!fH2D.count(h_name)) return;
+
+ auto h = fH2D[h_name].get();
+
+ for (int adc = 1; adc <= 31; adc++) {
+ TH1D* h_x = h->ProjectionX(Form("%s_dt", h_name.c_str()), adc, adc);
+
+ FitStsTime offset_fit(h_x);
+ TCanvas* drawing = offset_fit.Draw();
+ if (drawing != nullptr) {
+ std::string fit_folder = Form("0x%x/%02d", module, asic_idx);
+ TDirectory* dir = (TDirectory*) fReportFit->Get(fit_folder.c_str());
+ if (!dir) {
+ dir = fReportFit->mkdir(fit_folder.c_str());
+ }
+ dir->cd();
+ drawing->Write(Form("0x%x_%02d_%02d", module, asic_idx, adc));
+ }
+
+ if (offset_fit.use_fit_result) {
+ time_offset[n_pts] = offset_fit.fit_ptr->Parameter(1);
+ time_offset_err[n_pts] = offset_fit.fit_ptr->Parameter(2);
+ }
+ else {
+ /* If the fitting for individual ADC of the ASIC was unsuccessful
+ * the time offset is taken as the one estimated for the whole module
+ */
+ time_offset[n_pts] = module_offset;
+ time_offset_err[n_pts] = h_x->GetRMS();
+ }
+ adc_channel[n_pts] = adc;
+ adc_channel_err[n_pts] = 0;
+
+ // Update time walk map
+ fWalkMapRaw[module][asic_idx][adc] += time_offset[n_pts];
+
+ n_pts++;
+
+ } // ---- end ADC loop ----
+
+ fG1D[h_name] = std::make_unique<TGraphErrors>(n_pts, adc_channel, time_offset, adc_channel_err, time_offset_err);
+ fG1D[h_name]->SetTitle(Form("%s : StsDigi_0x%x time off-set", ToString(fRefSystem).c_str(), module));
+ fG1D[h_name]->GetXaxis()->SetTitle("Charge_{StsDigi} [ADC]");
+ fG1D[h_name]->GetYaxis()->SetTitle(Form("t_{%s} - t_{StsDigi} [ns]", ToString(fRefSystem).c_str()));
+ fG1D[h_name]->SetMinimum(fTimeWindowMin);
+ fG1D[h_name]->SetMaximum(fTimeWindowMax);
+
+ } // ---- end module loop ----
+ } // ---- end ASICs loop ----
+}
+
+double CbmStsTimeCal::FindGlobalOffset()
+{
+ double setup_offset = 0;
+ double min_sum_sigma = 999;
+ for (int32_t module : fAddressBook) { // module loop
+ double module_avg_sigma = 0;
+ double module_avg_offset = 0;
+ int n_of_values = 0;
+
+ for (unsigned int asic_idx = 0; asic_idx < 16; asic_idx++) { // asic loop
+ std::string h_name = Form("0x%x_%02d", module, asic_idx);
+ auto g = fG1D[h_name].get();
+
+ int n_of_pto = g->GetN();
+ int lower_bound = n_of_pto > 5 ? n_of_pto - 5 : 0;
+
+ // Calculate the average value of the sigma for the last ADC
+ for (int pto = n_of_pto; pto > lower_bound - 5; pto--) {
+ double offset = g->GetPointY(pto);
+ double sigma = g->GetErrorY(pto);
+ if (sigma) {
+ module_avg_sigma += sigma;
+ module_avg_offset += offset;
+ n_of_values++;
+ }
+ }
+ } // end asic loop
+
+ module_avg_sigma = n_of_values ? module_avg_sigma / n_of_values : -1;
+ module_avg_offset = n_of_values ? module_avg_offset / n_of_values : -1;
+
+ if (module_avg_sigma && module_avg_sigma < min_sum_sigma) {
+ min_sum_sigma = module_avg_sigma;
+ setup_offset = module_avg_offset;
+ }
+ } // end module loop
+
+ return setup_offset;
+}
+
+void CbmStsTimeCal::WriteYAML()
+{
+ std::string o_yaml_file = fOutputPath + "/StsTimeCalibration.yaml";
+ LOG(info) << "Writing parameter file to:" << o_yaml_file;
+
+ YAML::Emitter walk_map;
+ walk_map << YAML::BeginMap;
+ walk_map << YAML::Key << "timeOffset";
+ walk_map << YAML::Value << (int) fGlobalTimeOffset;
+ walk_map << YAML::Key << "WalkMap";
+ walk_map << YAML::Value << YAML::BeginMap;
+ for (const auto& [module, asics] : fWalkMapRaw) {
+ walk_map << YAML::Key << fmt::format("0x{:x}", module);
+ walk_map << YAML::Value << YAML::BeginSeq;
+ for (const auto& [asic, pars] : asics) {
+ // Begin a flow sequence
+ walk_map << YAML::Flow << YAML::BeginSeq;
+ for (const auto& value : pars) {
+ walk_map << fmt::format("{:+.3f}", value); // Add '+' for positive numbers
+ }
+ // End the flow sequence
+ walk_map << YAML::EndSeq;
+ }
+ walk_map << YAML::EndSeq;
+ }
+ walk_map << YAML::EndMap << YAML::EndMap;
+ assert(walk_map.good());
+
+ std::ofstream walk_map_out(o_yaml_file);
+ assert(walk_map_out.is_open());
+
+ walk_map_out << walk_map.c_str();
+ walk_map_out.close();
+}
+
+void CbmStsTimeCal::Finish()
+{
+ CheckTimeWalk();
+
+ // resize ADC offset vector to 31 channels from 0 [0-30`]
+ for (auto& [module, asics] : fWalkMapRaw) { // Loop over modules
+ for (auto& [asic_idx, asic_par] : asics) { // Loop over ASICs
+ for (int adc = 1; adc <= 31; adc++) { // loop over ADC channels 31 channels [1 - 31]
+ asic_par[adc - 1] = asic_par[adc];
+ }
+ asic_par.pop_back();
+ }
+ }
+
+ WriteYAML();
+
+ DrawResults();
+ SaveToFile();
+
+
+ fH1D.clear();
+ fH2D.clear();
+}
+
+
+void CbmStsTimeCal::DrawResults()
+{
+ std::unique_ptr<TCanvas> adc_vs_dt_1d =
+ std::make_unique<TCanvas>("adc_vs_dt_1d", "adc_vs_dt_1d", 10, 10, 4 * 1024, 4 * 1024);
+ adc_vs_dt_1d->Divide(4, 4);
+
+ std::unique_ptr<TCanvas> adc_vs_dt_2d =
+ std::make_unique<TCanvas>("adc_vs_dt_2d", "adc_vs_dt_2d", 10, 10, 4 * 1024, 4 * 1024);
+ adc_vs_dt_2d->Divide(4, 4);
+
+ std::unique_ptr<TCanvas> dt_vs_chn = std::make_unique<TCanvas>("dt_vs_chn", "dt_vs_chn", 10, 10, 1024, 1024);
+
+ fReportFile->cd();
+ for (int32_t sensor : fAddressBook) {
+ for (int asic_idx = 0; asic_idx < 16; asic_idx++) {
+ adc_vs_dt_1d->cd(asic_idx + 1);
+ if (fG1D[Form("0x%x_%02d", sensor, asic_idx)] != nullptr) {
+ fG1D[Form("0x%x_%02d", sensor, asic_idx)]->Draw("APL");
+ }
+
+ adc_vs_dt_2d->cd(asic_idx + 1);
+ if (fH2D[Form("0x%x_%02d", sensor, asic_idx)] != nullptr) {
+ fH2D[Form("0x%x_%02d", sensor, asic_idx)]->Draw("colz");
+ }
+ }
+
+ dt_vs_chn->cd();
+ if (fH2D[Form("0x%x_dt_vs_channel", sensor)] != nullptr) {
+ fH2D[Form("0x%x_dt_vs_channel", sensor)]->Draw("colz");
+ }
+
+ adc_vs_dt_1d->Write(Form("adc_vs_dt_1d_0x%x", sensor));
+ adc_vs_dt_2d->Write(Form("adc_vs_dt_2d_0x%x", sensor));
+ dt_vs_chn->Write(Form("dt_vs_chn_0x%x", sensor));
+ }
+}
diff --git a/analysis/detectors/sts/CbmStsTimeCal.h b/analysis/detectors/sts/CbmStsTimeCal.h
new file mode 100644
index 0000000000000000000000000000000000000000..a308ec8dfa6c3972cfc76b09083b60e98ed9174d
--- /dev/null
+++ b/analysis/detectors/sts/CbmStsTimeCal.h
@@ -0,0 +1,170 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+#ifndef CBMSTSTIMECAL_H
+#define CBMSTSTIMECAL_H
+
+#include "CbmDefs.h"
+#include "CbmDigiManager.h"
+#include "CbmStsAnaBase.h"
+#include "CbmStsUtils.h"
+#include "sts/WalkMap.h"
+#include "yaml/Yaml.h"
+
+#include <FairTask.h>
+
+#include <TClonesArray.h>
+
+#include <cstring>
+#include <map>
+#include <unordered_set>
+
+/** \class CbmStsTimeCal
+ * \brief Task for time calibration of STS hits.
+ *
+ * This class inherits from FairTask and CbmStsAnaBase.
+ * It provides functionality for calibrating the time of STS hits.
+ */
+
+/**
+ * @class CbmStsTimeCal
+ * @brief STS Digi-level time calibration task for ASIC and ADC synchronization.
+ *
+ * This class derives from FairTask and CbmStsAnaBase, and performs
+ * time calibration of the STS (Silicon Tracking System) by determining
+ * precise time offsets for each ASIC and ADC. The calibration removes
+ * time walk effects to improve timing accuracy.
+ *
+ * Key features:
+ *
+ * - **Time Offset Calculation**:
+ * Individual time offset values are found from the time difference distribution t_StsDigi - t_RefDigi (whole combinatoric inside search window).
+ * 1. Pre-determination of the time offset is obtained as the mean value of the time correlatiion peak for a given module
+ * 2. Further refining is attempted by fitting the time difference distribution with a N(dt,mu, sigma) + P_2(dt, p0,p1,p2).
+ * 3. ASIC ADC granularity is attemped. If failed, module offset is used.
+ *
+ * - **Reference Detector Selection**:
+ * A reference detector system can be specified.
+ * By default, the BMon (Beam Monitor) system is used.
+ *
+ * - **Reporting and Diagnostics**:
+ * Generates a set of reporting histograms that illustrate the timing
+ * distribution and calibration results.
+ *
+ * - **Organized Output by Run ID**:
+ * When a run ID is provided, all output files (e.g., parameters,
+ * histograms, reports) are stored in a dedicated folder named after
+ * the run ID, ensuring organized data handling.
+ *
+ * - **Output File Generation**:
+ * The task generates a YAML parameter file containing the individual time offset
+ * for each ASIC ADC.
+ *
+ * - **Using pre-existing time calibration**:
+ * Pre-existing YAML parameter files should be provided IF THE INPUT DATA WAS UNPACKED using such file.
+ *
+ * - **Time windows**:
+ * The time window for the calibration can be set using the lower and upper limits.
+ * - Very large windows can slow down the code due to high combinatorics.
+ * - To smmal windows can lead to biased result. The search windows should not be smaller than 25 ns.
+ *
+ * This task is essential for achieving precise timing synchronization
+ * across all readout electronics in the STS subsystem.
+ */
+class CbmStsTimeCal : public FairTask, public CbmStsAnaBase {
+ public:
+ CbmStsTimeCal() = default;
+ ~CbmStsTimeCal() = default;
+
+ /** \brief Parameterized constructor
+ * @param system reference detector for time
+ * @param lower_lim lower boundary for time difference between @param system digis and STS digis (default: -60).
+ * @param upper_lim lower boundary for time difference between @param system digis and STS digis (default: +60).
+ */
+ CbmStsTimeCal(ECbmModuleId system, double lower_lim = -60, double upper_lim = +60);
+
+ /** \brief Parameterized constructor
+ * @param system reference detector for time
+ * @param lower_lim lower boundary for time difference between @param system digis and STS digis (default: -60).
+ * @param upper_lim lower boundary for time difference between @param system digis and STS digis (default: +60).
+ */
+ CbmStsTimeCal(int run_id, ECbmModuleId system, double lower_lim, double upper_lim);
+
+ /**
+ * @brief Set parameter file used during unpacking
+ * @param file_name full/relative path to the parameters file
+ */
+ void SetWalkFile(std::string par_file);
+
+ InitStatus Init();
+ void Exec(Option_t*);
+ void Finish();
+
+
+ private:
+ double fGlobalTimeOffset{0};
+ double fTimeWindowMin{-60};
+ double fTimeWindowMax{+60};
+ const double kStsClock{3.125};
+ std::string fParFile{""};
+ std::string fOutputPath{""};
+
+ std::map<int /* Module */, std::map<int /* ASIC */, std::vector<double> /* ADC offset */>> fWalkMapRaw;
+ cbm::algo::sts::WalkMap fWalkMap;
+
+ std::unique_ptr<TFile> fReportOffset;
+ std::unique_ptr<TFile> fReportFit;
+
+ ECbmModuleId fRefSystem{ECbmModuleId::kBmon};
+
+ CbmDigiManager* fDigiManager{nullptr};
+
+ TClonesArray* fCbmEvtArray{nullptr};
+
+ /**
+ * @brief Book histograms.
+ * @param address The module address for which histograms will be booked.
+ */
+ void BookHistograms(int32_t address);
+
+ /** \brief Extract Time Walk parameters */
+ void CheckTimeWalk();
+
+ /** \brief Extract Time general offset for a given module */
+ double FindModuleOffset(int32_t);
+
+ /** \brief Find a common offset for STS setup.
+ * It is taken as the average of the time offset for ADC > 25
+ * for the setup module with the narrower dt distribution
+ */
+ double FindGlobalOffset();
+
+ /** \brief Check timing for a specific digi type
+ * It fills StsDigi time difference respect to fRefSystem Digis histograms
+ */
+ template<class Digi>
+ void CheckTiming();
+
+ /**
+ * @brief Initialized the Time Walk map
+ * @param address The module address for which the initialization is done
+ */
+ void InitTimeWalkMap(int32_t address);
+
+ /**
+ * @brief Load the time calibration parameters from user define file
+ */
+ void LoadWalkFromFile();
+
+ /**
+ * @brief Write parameters file in YAML format
+ * @param o_path output folder where to save the file
+ */
+ void WriteYAML();
+
+ void DrawResults();
+
+ ClassDef(CbmStsTimeCal, 1);
+};
+#endif
diff --git a/analysis/detectors/sts/CbmStsUtils.cxx b/analysis/detectors/sts/CbmStsUtils.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..c3b0c0c1bbca77d3b60f6158a1b123d1a98a40c5
--- /dev/null
+++ b/analysis/detectors/sts/CbmStsUtils.cxx
@@ -0,0 +1,143 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+#include "CbmStsUtils.h"
+
+#include "CbmStsAddress.h"
+#include "CbmStsCluster.h"
+#include "CbmStsParSetModule.h"
+
+#include <cassert>
+
+int32_t cbm_sts_utils::GetHitCluSizeF(CbmStsHit* hit, TClonesArray* sts_clu_array)
+{
+ if (hit == nullptr || sts_clu_array == nullptr) {
+ return -1;
+ }
+ return ((CbmStsCluster*) sts_clu_array->UncheckedAt(hit->GetFrontClusterId()))->GetSize();
+}
+int32_t cbm_sts_utils::GetHitCluSizeB(CbmStsHit* hit, TClonesArray* sts_clu_array)
+{
+ if (hit == nullptr || sts_clu_array == nullptr) {
+ return -1;
+ }
+ return ((CbmStsCluster*) sts_clu_array->UncheckedAt(hit->GetBackClusterId()))->GetSize();
+}
+/* Get StsHit charge as average of front and back cluster charges */
+double cbm_sts_utils::GetHitChargeF(CbmStsHit* hit, TClonesArray* sts_clu_array)
+{
+ if (hit == nullptr || sts_clu_array == nullptr) {
+ return -1;
+ }
+ return ((CbmStsCluster*) sts_clu_array->UncheckedAt(hit->GetFrontClusterId()))->GetCharge();
+}
+
+double cbm_sts_utils::GetHitChargeB(CbmStsHit* hit, TClonesArray* sts_clu_array)
+{
+ if (hit == nullptr || sts_clu_array == nullptr) {
+ return -1;
+ }
+ return ((CbmStsCluster*) sts_clu_array->UncheckedAt(hit->GetBackClusterId()))->GetCharge();
+}
+
+/* Get StsHit cluster time */
+double cbm_sts_utils::GetHitTimeF(CbmStsHit* hit, TClonesArray* sts_clu_array)
+{
+ if (hit == nullptr || sts_clu_array == nullptr) {
+ return -1;
+ }
+ return ((CbmStsCluster*) sts_clu_array->UncheckedAt(hit->GetFrontClusterId()))->GetTime();
+}
+
+double cbm_sts_utils::GetHitTimeB(CbmStsHit* hit, TClonesArray* sts_clu_array)
+{
+ if (hit == nullptr || sts_clu_array == nullptr) {
+ return -1;
+ }
+ return ((CbmStsCluster*) sts_clu_array->UncheckedAt(hit->GetBackClusterId()))->GetTime();
+}
+/* ---- ------------- ---- */
+
+double cbm_sts_utils::GetHitCharge(CbmStsHit* hit, TClonesArray* sts_clu_array)
+{
+ if (hit == nullptr || sts_clu_array == nullptr) {
+ return -1;
+ }
+ double cluster_charge_f = GetHitChargeF(hit, sts_clu_array);
+ double cluster_charge_b = GetHitChargeB(hit, sts_clu_array);
+ return 0.5 * (cluster_charge_f + cluster_charge_b);
+}
+
+/* Get charge asymmetry */
+double cbm_sts_utils::GetHitChargeAsy(CbmStsHit* hit, TClonesArray* sts_clu_array)
+{
+ if (hit == nullptr || sts_clu_array == nullptr) {
+ return -1;
+ }
+ double cluster_charge_f = GetHitChargeF(hit, sts_clu_array);
+ double cluster_charge_b = GetHitChargeB(hit, sts_clu_array);
+ return (cluster_charge_f - cluster_charge_b) / (cluster_charge_f + cluster_charge_b);
+}
+
+/** \brief Get the charge asymmetry of a hit **/
+std::pair<cbm_sts_utils::HBinning, cbm_sts_utils::HBinning>
+cbm_sts_utils::ChargeBinning(const CbmStsParModule& par_module, const uint32_t max_clu_size)
+{
+ assert(("Maximum cluster size must be greater than 0", max_clu_size > 0));
+
+ auto par_asic = par_module.GetAsicParams();
+
+ double n_side_q_thr = par_asic[0].GetThreshold();
+ double n_side_q_dyn = par_asic[0].GetDynRange();
+ double p_side_q_thr = par_asic[8].GetThreshold();
+ double p_side_q_dyn = par_asic[8].GetDynRange();
+
+ for (int asic_idx = 0; asic_idx < 16; asic_idx++) {
+ auto asic = par_asic[asic_idx];
+ if (asic_idx < 7) { // n-side
+ n_side_q_thr = std::min(n_side_q_thr, asic.GetThreshold());
+ n_side_q_dyn = std::max(n_side_q_dyn, asic.GetDynRange());
+ }
+ else { // p-side
+ p_side_q_thr = std::min(p_side_q_thr, asic.GetThreshold());
+ p_side_q_dyn = std::max(p_side_q_dyn, asic.GetDynRange());
+ }
+ }
+
+ double n_side_q_bin = n_side_q_dyn / 31;
+ double p_side_q_bin = p_side_q_dyn / 31;
+
+ double n_side_q_min = n_side_q_thr;
+ double p_side_q_min = p_side_q_thr;
+
+ double n_side_q_max = (n_side_q_dyn * (31.5 / 31) + 0.5 * n_side_q_bin) * max_clu_size + n_side_q_min;
+ double p_side_q_max = (p_side_q_dyn * (31.5 / 31) + 0.5 * p_side_q_bin) * max_clu_size + p_side_q_min;
+
+ uint32_t n_side_n_bin = 32 * max_clu_size;
+ uint32_t p_side_n_bin = 32 * max_clu_size;
+
+ return std::make_pair(cbm_sts_utils::HBinning{n_side_n_bin, n_side_q_min, n_side_q_max},
+ cbm_sts_utils::HBinning{p_side_n_bin, p_side_q_min, p_side_q_max});
+}
+
+
+std::set<int> cbm_sts_utils::GetUnits(const std::vector<int32_t> addresses)
+{
+ std::set<int> units;
+ for (int32_t address : addresses) {
+ units.insert(CbmStsAddress::GetElementId(address, kStsUnit));
+ }
+ return units;
+}
+
+std::vector<int32_t> cbm_sts_utils::GetUnitModules(const std::vector<int32_t> addresses, const uint32_t unit)
+{
+ std::vector<int32_t> unit_modules;
+ for (int32_t address : addresses) {
+ if (CbmStsAddress::GetElementId(address, kStsUnit) == unit) {
+ unit_modules.push_back(address);
+ }
+ }
+ return unit_modules;
+}
diff --git a/analysis/detectors/sts/CbmStsUtils.h b/analysis/detectors/sts/CbmStsUtils.h
new file mode 100644
index 0000000000000000000000000000000000000000..f6c4c98e16862e9fb460a7b0e411e9652fd0dbd6
--- /dev/null
+++ b/analysis/detectors/sts/CbmStsUtils.h
@@ -0,0 +1,118 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+#ifndef CBMSTSUTILS_H
+#define CBMSTSUTILS_H
+
+#include "CbmStsHit.h"
+#include "TClonesArray.h"
+
+#include <set>
+
+class CbmStsParModule;
+
+/** \namespace cbm_sts_utils
+ * \brief Namespace containing utility functions for STS detector.
+ */
+namespace cbm_sts_utils
+{
+ static const double kStsClock = 3.125;
+ static const double kStsDx = 6. / 1024;
+ static const double kStsDy = kStsDx / std::tan(7.5 * TMath::Pi() / 180);
+ static const double kStsErrX = kStsDx / sqrt(12);
+ static const double kStsErrY = kStsDy / sqrt(12);
+
+ /** \brief Get the cluster size of a hit from the front side
+ * \param hit Pointer to the STS hit
+ * \param clusters Pointer to the TClonesArray of clusters
+ * \return Cluster size on the front side
+ */
+ int32_t GetHitCluSizeF(CbmStsHit* hit = nullptr, TClonesArray* clusters = nullptr);
+
+ /** \brief Get the cluster size of a hit from the back side
+ * \param hit Pointer to the STS hit
+ * \param clusters Pointer to the TClonesArray of clusters
+ * \return Cluster size on the back side
+ */
+ int32_t GetHitCluSizeB(CbmStsHit* hit = nullptr, TClonesArray* clusters = nullptr);
+
+ /** \brief Get the charge of a hit as the average of front and back cluster charges
+ * \param hit Pointer to the STS hit
+ * \param clusters Pointer to the TClonesArray of clusters
+ * \return Average charge of the hit
+ */
+ double GetHitCharge(CbmStsHit* hit = nullptr, TClonesArray* clusters = nullptr);
+
+ /** \brief Get the charge of a hit from the back side
+ * \param hit Pointer to the STS hit
+ * \param clusters Pointer to the TClonesArray of clusters
+ * \return Charge of the hit on the back side
+ */
+ double GetHitChargeB(CbmStsHit* hit = nullptr, TClonesArray* clusters = nullptr);
+
+ /** \brief Get the charge of a hit from the front side
+ * \param hit Pointer to the STS hit
+ * \param clusters Pointer to the TClonesArray of clusters
+ * \return Charge of the hit on the front side
+ */
+ double GetHitChargeF(CbmStsHit* hit = nullptr, TClonesArray* clusters = nullptr);
+
+ /** \brief Get the charge of a hit from the back side
+ * \param hit Pointer to the STS hit
+ * \param clusters Pointer to the TClonesArray of clusters
+ * \return Charge of the hit on the back side
+ */
+ double GetHitTimeB(CbmStsHit* hit = nullptr, TClonesArray* clusters = nullptr);
+
+ /** \brief Get the charge of a hit from the front side
+ * \param hit Pointer to the STS hit
+ * \param clusters Pointer to the TClonesArray of clusters
+ * \return Charge of the hit on the front side
+ */
+ double GetHitTimeF(CbmStsHit* hit = nullptr, TClonesArray* clusters = nullptr);
+
+ /** \brief Get the size asymmetry of a hit
+ * \param hit Pointer to the STS hit
+ * \param clusters Pointer to the TClonesArray of clusters
+ * \return Charge asymmetry of the hit
+ */
+ double GetHitSizeAsy(CbmStsHit* hit = nullptr, TClonesArray* clusters = nullptr);
+
+ /** \brief Get the charge asymmetry of a hit
+ * \param hit Pointer to the STS hit
+ * \param clusters Pointer to the TClonesArray of clusters
+ * \return Charge asymmetry of the hit
+ */
+ double GetHitChargeAsy(CbmStsHit* hit = nullptr, TClonesArray* clusters = nullptr);
+
+ /**
+ * \brief Structure to hold the binning for 1D histogram
+ **/
+ struct HBinning {
+ uint32_t n_of_bins{0};
+ double x_min{0};
+ double x_max{0};
+ };
+
+ /** \brief Generate the charge binning from module config obj
+ * \param par_module Obj containing the module configuration
+ * \param max_clu_size maximum cluster size to sample
+ * \return pair binning structure for n(p)-side
+ */
+ std::pair<HBinning, HBinning> ChargeBinning(const CbmStsParModule& par_module, const uint32_t max_clu_size = 1);
+
+ /** \brief Return the STS units from a list of addresses
+ * \param addresses Vector containing the addresses of the modules
+ * \return set of units
+ */
+ std::set<int> GetUnits(const std::vector<int32_t> addresses);
+
+ /** \brief From a list of module address, return those that belong to a selected unit
+ * \param addresses Vector containing the addresses of the modules
+ * \param unit Selected unit
+ * \return subset of modules that belong to the selected unit
+ */
+ std::vector<int32_t> GetUnitModules(const std::vector<int32_t> addresses, const uint32_t unit);
+}; // namespace cbm_sts_utils
+#endif
diff --git a/analysis/detectors/sts/tests/CMakeLists.txt b/analysis/detectors/sts/tests/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..49e69527ae7f3e0a1379adaf22f2cca5a38c7b7e
--- /dev/null
+++ b/analysis/detectors/sts/tests/CMakeLists.txt
@@ -0,0 +1,27 @@
+# CMakeList file for library libStsAna dedicated tests
+# Last update: Dario Ramirez, 23.04.2025
+
+Function(AddBasicTest name)
+ set(PVT_DEPS
+ Gtest
+ GtestMain
+ CbmData
+ CbmStsBase
+ CbmStsAna
+ )
+
+ add_executable(${name} ${name}.cxx)
+ target_include_directories(${name} PRIVATE ${INCLUDE_DIRECTORIES})
+ target_link_libraries(${name} PRIVATE ${PVT_DEPS})
+ Add_Test(
+ NAME ${name}
+ COMMAND ${name}
+ WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}
+ )
+EndFunction()
+
+set(INCLUDE_DIRECTORIES
+ ${CMAKE_SOURCE_DIR}/analysis/detectors/sts
+ )
+
+AddBasicTest(_GTestCbmCut)
diff --git a/analysis/detectors/sts/tests/_GTestCbmCut.cxx b/analysis/detectors/sts/tests/_GTestCbmCut.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..1db1ccf035e3e875ad8dfba553743dae077a20b3
--- /dev/null
+++ b/analysis/detectors/sts/tests/_GTestCbmCut.cxx
@@ -0,0 +1,62 @@
+/* Copyright (C) 2016-2020 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+#include "CbmCut.h"
+#include "gtest/gtest.h"
+
+TEST(_GTestCbmCut, DefConstructor)
+{
+ CbmCut<float> cut;
+ EXPECT_EQ(false, cut.GetMinState());
+ EXPECT_EQ(false, cut.GetMaxState());
+}
+
+TEST(_GTestCbmCut, ParConstructor)
+{
+ float min = -100;
+ float max = +100;
+ CbmCut<float> cut(min, max);
+ EXPECT_EQ(true, cut.GetMinState());
+ EXPECT_EQ(true, cut.GetMaxState());
+ EXPECT_EQ(min, cut.GetMin());
+ EXPECT_EQ(max, cut.GetMax());
+}
+
+TEST(_GTestCbmCut, CheckOnDefConstructor)
+{
+ CbmCut<float> cut;
+ for (float i = -1e+6; i < 1e+6; i++) {
+ EXPECT_EQ(true, cut.Check(i));
+ }
+}
+
+TEST(_GTestCbmCut, CheckCase1)
+{
+ float min = -100;
+ float max = +100;
+ CbmCut<float> cut(min, max);
+ for (float i = -1e+6; i < 1e+6; i++) {
+ if (i < min || i > max) {
+ EXPECT_EQ(false, cut.Check(i));
+ }
+ else {
+ EXPECT_EQ(true, cut.Check(i));
+ }
+ }
+}
+
+TEST(_GTestCbmCut, CheckCase2)
+{
+ float min = +100;
+ float max = -100;
+ CbmCut<float> cut(min, max);
+ for (float i = -1e+6; i < 1e+6; i++) {
+ if (i <= max || i >= min) {
+ EXPECT_EQ(true, cut.Check(i));
+ }
+ else {
+ EXPECT_EQ(false, cut.Check(i));
+ }
+ }
+}
diff --git a/core/detectors/sts/CbmStsParSetModule.cxx b/core/detectors/sts/CbmStsParSetModule.cxx
index b9439025c79e6912447ecb7383a087edaa0c6cb9..62b888a00d98379c01643e762ec469cff6c92c12 100644
--- a/core/detectors/sts/CbmStsParSetModule.cxx
+++ b/core/detectors/sts/CbmStsParSetModule.cxx
@@ -13,6 +13,7 @@
#include <Logger.h> // for LOG, Logger
#include <cassert> // for assert
+#include <fstream> // for reading parameters from ASCII file
#include <sstream> // for operator<<, basic_ostream, stringstream
#include <string> // for char_traits
@@ -64,6 +65,84 @@ Bool_t CbmStsParSetModule::getParams(FairParamList*)
}
// --------------------------------------------------------------------------
+// ----- Read parameters from ASCII file --------------------------------
+Bool_t CbmStsParSetModule::LoadParASCII(std::string file_name)
+{
+ std::ifstream input(file_name);
+ if (!input.is_open()) {
+ LOG(error) << Form("[STS Charge Calibration] Error opening file: %s", file_name.c_str());
+ ;
+ return kFALSE;
+ }
+
+ LOG(info) << Form("[STS Charge Calibration] Loading configuration from file: %s", file_name.c_str());
+
+ int32_t address, side, asic_idx, nb_channels, nb_adc;
+ double dyn_range, threshold, time_res, dead_time, noise, znr;
+
+ // Default ASIC values
+ std::unique_ptr<CbmStsParAsic> default_asic =
+ std::make_unique<CbmStsParAsic>(128, 31, 75000., 3000., 5., 800., 1000., 3.9789e-3);
+
+ // Default module configuration
+ std::unique_ptr<CbmStsParModule> default_module = std::make_unique<CbmStsParModule>(2048, 128);
+ default_module->SetAllAsics(*default_asic);
+
+ std::string line;
+ while (std::getline(input, line)) {
+
+ // Skip commented lines
+ if (line[0] == '#') {
+ continue;
+ }
+ std::istringstream str_stream(line);
+ str_stream >> std::hex >> address >> std::dec >> side >> asic_idx >> nb_channels >> nb_adc >> dyn_range >> threshold
+ >> time_res >> dead_time >> noise >> znr;
+
+ if (address == -1) {
+ LOG(info) << "[STS Charge Calibration] Using same configuration for all modules";
+ SetGlobalPar(*default_module);
+ break;
+ }
+
+ // Initialize map entry with default configuration
+ if (!fParams.count(address)) {
+ fParams[address] = *default_module;
+ }
+
+ // Custom ASIC loaded from file
+ std::unique_ptr<CbmStsParAsic> custom_asic =
+ std::make_unique<CbmStsParAsic>(nb_channels, nb_adc, dyn_range, threshold, time_res, dead_time, noise, znr);
+
+ if (asic_idx == -1) { // Same configuration for all ASIC in the given side
+ if (side == -1) { // Same configuration for both sides
+ LOG(info) << Form("[STS Charge Calibration] Using same configuration for all ASICs in module 0x%x", address);
+ fParams[address].SetAllAsics(*custom_asic);
+ }
+ else {
+ // Configure ASIC all for given side
+ LOG(info) << Form("[STS Charge Calibration] Using same configuration for all ASICs in module 0x%x, side %u",
+ address, side);
+ for (int idx = 0; idx < 8; idx++) {
+ fParams[address].SetAsic(idx + 8 * side, *custom_asic);
+ }
+ }
+ }
+ else {
+ if (side == -1) {
+ fParams[address].SetAsic(asic_idx, *custom_asic);
+ fParams[address].SetAsic(asic_idx + 8, *custom_asic);
+ }
+ else {
+ fParams[address].SetAsic(asic_idx + 8 * side, *custom_asic);
+ }
+ }
+ }
+
+
+ return kTRUE;
+}
+// --------------------------------------------------------------------------
// ----- Get condition parameters of a module ---------------------------
const CbmStsParModule& CbmStsParSetModule::GetParModule(UInt_t address)
diff --git a/core/detectors/sts/CbmStsParSetModule.h b/core/detectors/sts/CbmStsParSetModule.h
index 3f516fc4f4adadb24d882c1d29521141529a3c71..6302ee1d96376e02c17a9b9618e68ba9c10b129b 100644
--- a/core/detectors/sts/CbmStsParSetModule.h
+++ b/core/detectors/sts/CbmStsParSetModule.h
@@ -64,6 +64,11 @@ public:
**/
virtual Bool_t getParams(FairParamList* parList);
+ /** @brief Reading parameters from ASCII.
+ * * @param file_name path to the parameters file
+ **/
+ Bool_t LoadParASCII(std::string file_name);
+
/** @brief Get condition parameters of a sensor
** @param Module address
diff --git a/macro/sts/cbm_vertex.C b/macro/sts/cbm_vertex.C
new file mode 100755
index 0000000000000000000000000000000000000000..12b6baada433f388f5f3db027d90c5c4825cc315
--- /dev/null
+++ b/macro/sts/cbm_vertex.C
@@ -0,0 +1,101 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+void cbm_vertex(int n_of_ts = -1,
+ double dca_cut = 0.5, // cm
+ std::string geo_setup_tag = "mcbm_beam_2024_05_08_nickel", std::string raw_file = "",
+ std::string rec_file = "", std::string geo_file = "", std::string out_file = "cbm_vertex_dca.root",
+ bool bAli = true)
+{
+ // ----- Timer --------------------------------------------------------
+ TStopwatch timer;
+ timer.Start();
+
+ // --- Logger settings ----------------------------------------------------
+ TString logLevel = "DEBUG";
+ TString logVerbosity = "veryhigh";
+ logLevel = "INFO";
+ logVerbosity = "low";
+
+ FairLogger::GetLogger()->SetLogScreenLevel(logLevel.Data());
+ FairLogger::GetLogger()->SetLogVerbosityLevel(logVerbosity.Data());
+ // ------------------------------------------------------------------------
+
+
+ // ----- Environment --------------------------------------------------
+ std::string srcDir = gSystem->Getenv("VMCWORKDIR"); // top source directory
+ // ------------------------------------------------------------------------
+
+ std::unique_ptr<CbmStsAnalysis> run = std::make_unique<CbmStsAnalysis>();
+
+ FairFileSource* inputSource = new FairFileSource(rec_file.c_str());
+ inputSource->AddFriend(raw_file.c_str());
+
+ FairRootFileSink* outputSink = new FairRootFileSink(out_file.c_str());
+
+
+ if (bAli) {
+ TString alignmentMatrixFileName = srcDir + "/parameters/mcbm/AlignmentMatrices_" + geo_setup_tag + ".root";
+ if (alignmentMatrixFileName.Length() != 0) {
+ std::map<std::string, TGeoHMatrix>* matrices{nullptr};
+ TFile* misalignmentMatrixRootfile = new TFile(alignmentMatrixFileName, "READ");
+ if (misalignmentMatrixRootfile->IsOpen()) {
+ gDirectory->GetObject("MisalignMatrices", matrices);
+ misalignmentMatrixRootfile->Close();
+ }
+ else {
+ LOG(error) << "Could not open file " << alignmentMatrixFileName << "\n Exiting";
+ exit(1);
+ }
+
+ if (matrices) {
+ run->AddAlignmentMatrices(*matrices);
+ }
+ else {
+ LOG(error) << "Alignment required but no matrices found."
+ << "\n Exiting";
+ exit(1);
+ }
+ }
+ }
+
+ run->SetGeomFile(geo_file.c_str());
+ run->SetSource(inputSource);
+ run->SetSink(outputSink);
+
+
+ // ------------------------------------------------------------------------
+ // Configure Analysis filters
+ // ------------------------------------------------------------------------
+ CbmCutMap ana_filter;
+ ana_filter.AddCbmCut(CbmCutId::kGlobalTrackStsSize)->SetMin(2);
+ ana_filter.AddCbmCut(CbmCutId::kGlobalTrackTrdSize)->SetMin(1);
+ ana_filter.AddCbmCut(CbmCutId::kGlobalTrackTofSize)->SetMin(1);
+ // ana_filter.AddCbmCut(CbmCutId::kGlobalTrackPval)->SetMax(0.005);
+
+ // ------------------------------------------------------------------------
+ // Configure Vertex Finder
+ // ------------------------------------------------------------------------
+ std::shared_ptr<CbmDcaVertexFinder> vertex_finder = std::make_shared<CbmDcaVertexFinder>(dca_cut);
+
+ // ------------------------------------------------------------------------
+ // Configure Task
+ // ------------------------------------------------------------------------
+ CbmEventVertexDca* cbm_vtx_cda = new CbmEventVertexDca();
+ cbm_vtx_cda->SetVertexFinder(vertex_finder);
+ cbm_vtx_cda->SetCutMap(&ana_filter);
+ run->AddTask(cbm_vtx_cda);
+
+ run->Init();
+ run->Run(0, n_of_ts);
+
+ // ------------------------------------------------------------------------
+ timer.Stop();
+ Double_t rtime = timer.RealTime();
+ Double_t ctime = timer.CpuTime();
+ cout << endl << endl;
+ cout << "Macro finished successfully." << endl;
+ cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
+ cout << endl;
+}
diff --git a/macro/sts/sts_beam_spot.C b/macro/sts/sts_beam_spot.C
new file mode 100644
index 0000000000000000000000000000000000000000..c927fcadd7ab7587d021b702037e1e7f9bd8e073
--- /dev/null
+++ b/macro/sts/sts_beam_spot.C
@@ -0,0 +1,123 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+void sts_beam_spot(int n_of_ts = -1, std::string geo_setup_tag = "mcbm_beam_2024_05_08_nickel",
+ std::string raw_file = "", std::string rec_file = "", std::string geo_file = "",
+ std::string out_file = "sts_beam_spot.root", bool bAli = true)
+{
+ // ----- Timer --------------------------------------------------------
+ TStopwatch timer;
+ timer.Start();
+
+ // --- Logger settings ----------------------------------------------------
+ TString logLevel = "DEBUG";
+ TString logVerbosity = "veryhigh";
+ logLevel = "INFO";
+ logVerbosity = "low";
+
+ FairLogger::GetLogger()->SetLogScreenLevel(logLevel.Data());
+ FairLogger::GetLogger()->SetLogVerbosityLevel(logVerbosity.Data());
+ // ------------------------------------------------------------------------
+
+
+ // ----- Environment --------------------------------------------------
+ std::string srcDir = gSystem->Getenv("VMCWORKDIR"); // top source directory
+ // ------------------------------------------------------------------------
+
+ std::unique_ptr<CbmStsAnalysis> run = std::make_unique<CbmStsAnalysis>();
+
+ FairFileSource* inputSource = new FairFileSource(rec_file.c_str());
+ inputSource->AddFriend(raw_file.c_str());
+
+ FairRootFileSink* outputSink = new FairRootFileSink(out_file.c_str());
+
+
+ if (bAli) {
+ TString alignmentMatrixFileName = srcDir + "/parameters/mcbm/AlignmentMatrices_" + geo_setup_tag + ".root";
+ if (alignmentMatrixFileName.Length() != 0) {
+ std::map<std::string, TGeoHMatrix>* matrices{nullptr};
+ TFile* misalignmentMatrixRootfile = new TFile(alignmentMatrixFileName, "READ");
+ if (misalignmentMatrixRootfile->IsOpen()) {
+ gDirectory->GetObject("MisalignMatrices", matrices);
+ misalignmentMatrixRootfile->Close();
+ }
+ else {
+ LOG(error) << "Could not open file " << alignmentMatrixFileName << "\n Exiting";
+ exit(1);
+ }
+
+ if (matrices) {
+ run->AddAlignmentMatrices(*matrices);
+ }
+ else {
+ LOG(error) << "Alignment required but no matrices found."
+ << "\n Exiting";
+ exit(1);
+ }
+ }
+ }
+
+ run->SetGeomFile(geo_file.c_str());
+ run->SetSource(inputSource);
+ run->SetSink(outputSink);
+
+
+ // ------------------------------------------------------------------------
+ // Configure Analysis filters
+ // ------------------------------------------------------------------------
+ CbmCutMap ana_filter;
+ ana_filter.AddCbmCut(CbmCutId::kStsHitCharge)->SetMin(10000);
+ ana_filter.AddCbmCut(CbmCutId::kEventNofStsHit)->SetMin(2);
+
+ // ------------------------------------------------------------------------
+ // Configure mCBM target planes
+ // ------------------------------------------------------------------------
+ double l_k0 = 38.45;
+ double l_k1 = l_k0 + 0.86 + 11.8 + 1.3 + 0.7;
+ double l_b1 = l_k0 + 4.99 + 1.2 + 0.75;
+
+ double beam_rot_y = (25.) * TMath::DegToRad();
+ CbmTarget Ni;
+ Ni.SetPosition(0, 0, 0);
+ Ni.SetRotation(beam_rot_y);
+
+ CbmTarget T0;
+ T0.SetPosition(-20 * sin(beam_rot_y), 0, -20 * cos(beam_rot_y));
+ T0.SetRotation(beam_rot_y);
+
+ CbmTarget K0;
+ K0.SetPosition(-l_k0 * sin(beam_rot_y), 0, -l_k0 * cos(beam_rot_y));
+ K0.SetRotation(beam_rot_y);
+
+ CbmTarget B1;
+ B1.SetPosition(-l_b1 * sin(beam_rot_y), 0, -l_b1 * cos(beam_rot_y));
+ B1.SetRotation(beam_rot_y);
+
+ CbmTarget K1;
+ K1.SetPosition(-l_k1 * sin(beam_rot_y), 0, -l_k1 * cos(beam_rot_y));
+ K1.SetRotation(beam_rot_y);
+ // ------------------------------------------------------------------------
+
+ CbmStsRecoBeamSpot* sts_bs = new CbmStsRecoBeamSpot();
+ sts_bs->SetCutMap(&ana_filter);
+ sts_bs->AddTarget("Ni", &Ni);
+ sts_bs->AddTarget("T0", &T0);
+ sts_bs->AddTarget("K0", &K0);
+ sts_bs->AddTarget("B1", &B1);
+ sts_bs->AddTarget("K1", &K1);
+
+ run->AddTask(sts_bs);
+
+ run->Init();
+ run->Run(0, n_of_ts);
+
+ // ------------------------------------------------------------------------
+ timer.Stop();
+ Double_t rtime = timer.RealTime();
+ Double_t ctime = timer.CpuTime();
+ cout << endl << endl;
+ cout << "Macro finished successfully." << endl;
+ cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
+ cout << endl;
+}
diff --git a/macro/sts/sts_correlation.C b/macro/sts/sts_correlation.C
new file mode 100644
index 0000000000000000000000000000000000000000..ac5b37e11bfc4b91ec6e12d38289b6bc01805b1d
--- /dev/null
+++ b/macro/sts/sts_correlation.C
@@ -0,0 +1,91 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+void sts_correlation(int n_of_ts = -1, std::string geo_setup_tag = "mcbm_beam_2024_05_08_nickel",
+ std::string raw_file = "", std::string rec_file = "", std::string geo_file = "",
+ std::string out_file = "cbm_sts_correlation.root", bool bAli = true)
+{
+ // ----- Timer --------------------------------------------------------
+ TStopwatch timer;
+ timer.Start();
+
+ // --- Logger settings ----------------------------------------------------
+ TString logLevel = "DEBUG";
+ TString logVerbosity = "veryhigh";
+ logLevel = "INFO";
+ logVerbosity = "low";
+
+ FairLogger::GetLogger()->SetLogScreenLevel(logLevel.Data());
+ FairLogger::GetLogger()->SetLogVerbosityLevel(logVerbosity.Data());
+ // ------------------------------------------------------------------------
+
+
+ // ----- Environment --------------------------------------------------
+ std::string srcDir = gSystem->Getenv("VMCWORKDIR"); // top source directory
+ // ------------------------------------------------------------------------
+
+ std::unique_ptr<CbmStsAnalysis> run = std::make_unique<CbmStsAnalysis>();
+
+ FairFileSource* inputSource = new FairFileSource(rec_file.c_str());
+ inputSource->AddFriend(raw_file.c_str());
+
+ FairRootFileSink* outputSink = new FairRootFileSink(out_file.c_str());
+
+
+ if (bAli) {
+ TString alignmentMatrixFileName = srcDir + "/parameters/mcbm/AlignmentMatrices_" + geo_setup_tag + ".root";
+ if (alignmentMatrixFileName.Length() != 0) {
+ std::map<std::string, TGeoHMatrix>* matrices{nullptr};
+ TFile* misalignmentMatrixRootfile = new TFile(alignmentMatrixFileName, "READ");
+ if (misalignmentMatrixRootfile->IsOpen()) {
+ gDirectory->GetObject("MisalignMatrices", matrices);
+ misalignmentMatrixRootfile->Close();
+ }
+ else {
+ LOG(error) << "Could not open file " << alignmentMatrixFileName << "\n Exiting";
+ exit(1);
+ }
+
+ if (matrices) {
+ run->AddAlignmentMatrices(*matrices);
+ }
+ else {
+ LOG(error) << "Alignment required but no matrices found."
+ << "\n Exiting";
+ exit(1);
+ }
+ }
+ }
+
+ run->SetGeomFile(geo_file.c_str());
+ run->SetSource(inputSource);
+ run->SetSink(outputSink);
+
+
+ // ------------------------------------------------------------------------
+ // Configure Analysis filters
+ // ------------------------------------------------------------------------
+ CbmCutMap ana_filter;
+ ana_filter.AddCbmCut(CbmCutId::kStsHitCharge)->SetMin(10000);
+ ana_filter.AddCbmCut(CbmCutId::kEventNofStsHit)->SetMin(2);
+
+ // ------------------------------------------------------------------------
+ // Configure Correlation
+ // ------------------------------------------------------------------------
+ CbmStsCorrelation* sts_corr = new CbmStsCorrelation();
+ sts_corr->SetCutMap(&ana_filter);
+ run->AddTask(sts_corr);
+
+ run->Init();
+ run->Run(0, n_of_ts);
+
+ // ------------------------------------------------------------------------
+ timer.Stop();
+ Double_t rtime = timer.RealTime();
+ Double_t ctime = timer.CpuTime();
+ cout << endl << endl;
+ cout << "Macro finished successfully." << endl;
+ cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
+ cout << endl;
+}
diff --git a/macro/sts/sts_efficiency.C b/macro/sts/sts_efficiency.C
new file mode 100755
index 0000000000000000000000000000000000000000..091c43dbc4c9de3d30542c422a4dcb61f49900fd
--- /dev/null
+++ b/macro/sts/sts_efficiency.C
@@ -0,0 +1,107 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+void sts_efficiency(int n_of_ts = -1, std::string geo_setup_tag = "mcbm_beam_2024_05_08_nickel",
+ std::string raw_file = "", std::string rec_file = "", std::string geo_file = "",
+ std::string out_file = "cbm_sts_efficiency.root", bool bAli = true)
+{
+ // ----- Timer --------------------------------------------------------
+ TStopwatch timer;
+ timer.Start();
+
+ // --- Logger settings ----------------------------------------------------
+ TString logLevel = "DEBUG";
+ TString logVerbosity = "veryhigh";
+ logLevel = "INFO";
+ logVerbosity = "low";
+
+ FairLogger::GetLogger()->SetLogScreenLevel(logLevel.Data());
+ FairLogger::GetLogger()->SetLogVerbosityLevel(logVerbosity.Data());
+ // ------------------------------------------------------------------------
+
+
+ // ----- Environment --------------------------------------------------
+ std::string srcDir = gSystem->Getenv("VMCWORKDIR"); // top source directory
+ // ------------------------------------------------------------------------
+
+ std::unique_ptr<CbmStsAnalysis> run = std::make_unique<CbmStsAnalysis>();
+
+ FairFileSource* inputSource = new FairFileSource(rec_file.c_str());
+ inputSource->AddFriend(raw_file.c_str());
+
+ FairRootFileSink* outputSink = new FairRootFileSink(out_file.c_str());
+
+
+ if (bAli) {
+ TString alignmentMatrixFileName = srcDir + "/parameters/mcbm/AlignmentMatrices_" + geo_setup_tag + ".root";
+ if (alignmentMatrixFileName.Length() != 0) {
+ std::map<std::string, TGeoHMatrix>* matrices{nullptr};
+ TFile* misalignmentMatrixRootfile = new TFile(alignmentMatrixFileName, "READ");
+ if (misalignmentMatrixRootfile->IsOpen()) {
+ gDirectory->GetObject("MisalignMatrices", matrices);
+ misalignmentMatrixRootfile->Close();
+ }
+ else {
+ LOG(error) << "Could not open file " << alignmentMatrixFileName << "\n Exiting";
+ exit(1);
+ }
+
+ if (matrices) {
+ run->AddAlignmentMatrices(*matrices);
+ }
+ else {
+ LOG(error) << "Alignment required but no matrices found."
+ << "\n Exiting";
+ exit(1);
+ }
+ }
+ }
+
+ run->SetGeomFile(geo_file.c_str());
+ run->SetSource(inputSource);
+ run->SetSink(outputSink);
+
+
+ // ------------------------------------------------------------------------
+ // Configure Analysis filters
+ // ------------------------------------------------------------------------
+ CbmCutMap ana_filter;
+ ana_filter.AddCbmCut(CbmCutId::kStsHitCharge)->SetMin(10000);
+ ana_filter.AddCbmCut(CbmCutId::kEventNofStsHit)->SetMin(2);
+ ana_filter.AddCbmCut(CbmCutId::kGlobalTrackStsSize)->SetMin(2);
+ ana_filter.AddCbmCut(CbmCutId::kGlobalTrackTrdSize)->SetMin(1);
+ ana_filter.AddCbmCut(CbmCutId::kGlobalTrackTofSize)->SetMin(1);
+ // ana_filter.AddCbmCut(CbmCutId::kGlobalTrackPval)->SetMax(0.005);
+
+ // ------------------------------------------------------------------------
+ // Configure Vertex Finder
+ // ------------------------------------------------------------------------
+ std::shared_ptr<CbmDcaVertexFinder> vertex_finder = std::make_shared<CbmDcaVertexFinder>(0.2);
+
+ // ------------------------------------------------------------------------
+ // Configure Efficiency
+ // ------------------------------------------------------------------------
+ CbmStsEfficiency* sts_eff1 = new CbmStsEfficiency(1, // test_layer_(layer_idx)
+ 0.25, // max_dis_trg_vtx_(cut_trg_vtx)
+ 0.25, // max_dca_trk_vtx_(cut_trk_vtx)
+ 0.005 // default_residual(def_res)
+ );
+
+ sts_eff1->SetCutMap(&ana_filter);
+ sts_eff1->SetVertexFinder(vertex_finder.get());
+ // sts_eff1->SetResidual("Ref_Sts10_residuals.info");
+ run->AddTask(sts_eff1);
+
+ run->Init();
+ run->Run(0, n_of_ts);
+
+ // ------------------------------------------------------------------------
+ timer.Stop();
+ Double_t rtime = timer.RealTime();
+ Double_t ctime = timer.CpuTime();
+ cout << endl << endl;
+ cout << "Macro finished successfully." << endl;
+ cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
+ cout << endl;
+}
diff --git a/macro/sts/sts_hit_ana.C b/macro/sts/sts_hit_ana.C
new file mode 100644
index 0000000000000000000000000000000000000000..327d69b2ba271f58e3f8212a876e9ba1a44c508f
--- /dev/null
+++ b/macro/sts/sts_hit_ana.C
@@ -0,0 +1,97 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+void sts_hit_ana(int n_of_ts = -1, std::string geo_setup_tag = "mcbm_beam_2024_05_08_nickel", std::string raw_file = "",
+ std::string rec_file = "", std::string geo_file = "", std::string out_file = "cbm_sts_hit_ana.root",
+ bool bAli = true)
+{
+ // ----- Timer --------------------------------------------------------
+ TStopwatch timer;
+ timer.Start();
+
+ // --- Logger settings ----------------------------------------------------
+ TString logLevel = "DEBUG";
+ TString logVerbosity = "veryhigh";
+ logLevel = "INFO";
+ logVerbosity = "low";
+
+ FairLogger::GetLogger()->SetLogScreenLevel(logLevel.Data());
+ FairLogger::GetLogger()->SetLogVerbosityLevel(logVerbosity.Data());
+ // ------------------------------------------------------------------------
+
+
+ // ----- Environment --------------------------------------------------
+ std::string srcDir = gSystem->Getenv("VMCWORKDIR"); // top source directory
+ // ------------------------------------------------------------------------
+
+ std::unique_ptr<CbmStsAnalysis> run = std::make_unique<CbmStsAnalysis>();
+
+ FairFileSource* inputSource = new FairFileSource(rec_file.c_str());
+ inputSource->AddFriend(raw_file.c_str());
+
+ FairRootFileSink* outputSink = new FairRootFileSink(out_file.c_str());
+
+
+ if (bAli) {
+ TString alignmentMatrixFileName = srcDir + "/parameters/mcbm/AlignmentMatrices_" + geo_setup_tag + ".root";
+ if (alignmentMatrixFileName.Length() != 0) {
+ std::map<std::string, TGeoHMatrix>* matrices{nullptr};
+ TFile* misalignmentMatrixRootfile = new TFile(alignmentMatrixFileName, "READ");
+ if (misalignmentMatrixRootfile->IsOpen()) {
+ gDirectory->GetObject("MisalignMatrices", matrices);
+ misalignmentMatrixRootfile->Close();
+ }
+ else {
+ LOG(error) << "Could not open file " << alignmentMatrixFileName << "\n Exiting";
+ exit(1);
+ }
+
+ if (matrices) {
+ run->AddAlignmentMatrices(*matrices);
+ }
+ else {
+ LOG(error) << "Alignment required but no matrices found."
+ << "\n Exiting";
+ exit(1);
+ }
+ }
+ }
+
+ run->SetGeomFile(geo_file.c_str());
+ run->SetSource(inputSource);
+ run->SetSink(outputSink);
+
+
+ // ------------------------------------------------------------------------
+ // Configure Analysis filters
+ // ------------------------------------------------------------------------
+ CbmCutMap ana_filter;
+ // Filter for all STS hits
+ ana_filter.AddCbmCut(CbmCutId::kStsHitQasym)->SetRange(-0.8, +0.8);
+ ana_filter.AddCbmCut(CbmCutId::kStsHitCharge)->SetMin(10000);
+
+ // Filter for track to extract STS hits tracks
+ ana_filter.AddCbmCut(CbmCutId::kGlobalTrackStsSize)->SetMin(2);
+ ana_filter.AddCbmCut(CbmCutId::kGlobalTrackTrdSize)->SetMin(1);
+ ana_filter.AddCbmCut(CbmCutId::kGlobalTrackTofSize)->SetMin(1);
+
+ // ------------------------------------------------------------------------
+ // Configure Charge Analysis
+ // ------------------------------------------------------------------------
+ CbmStsHitAna* sts_hit = new CbmStsHitAna(Form("%s/parameters/sts/sts_charge_cal_v24a.par", srcDir.c_str()));
+ sts_hit->SetCutMap(&ana_filter);
+ run->AddTask(sts_hit);
+
+ run->Init();
+ run->Run(0, n_of_ts);
+
+ // ------------------------------------------------------------------------
+ timer.Stop();
+ Double_t rtime = timer.RealTime();
+ Double_t ctime = timer.CpuTime();
+ cout << endl << endl;
+ cout << "Macro finished successfully." << endl;
+ cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
+ cout << endl;
+}
diff --git a/macro/sts/sts_resolution.C b/macro/sts/sts_resolution.C
new file mode 100755
index 0000000000000000000000000000000000000000..75404c692627d912190878cc70de113fa94c357d
--- /dev/null
+++ b/macro/sts/sts_resolution.C
@@ -0,0 +1,98 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+void sts_resolution(int n_of_ts = -1, std::string geo_setup_tag = "mcbm_beam_2024_05_08_nickel",
+ std::string raw_file = "", std::string rec_file = "", std::string geo_file = "",
+ std::string out_file = "cbm_sts_correlation.root", bool bAli = true)
+{
+ // ----- Timer --------------------------------------------------------
+ TStopwatch timer;
+ timer.Start();
+
+ // --- Logger settings ----------------------------------------------------
+ TString logLevel = "DEBUG";
+ TString logVerbosity = "veryhigh";
+ logLevel = "INFO";
+ logVerbosity = "low";
+
+ FairLogger::GetLogger()->SetLogScreenLevel(logLevel.Data());
+ FairLogger::GetLogger()->SetLogVerbosityLevel(logVerbosity.Data());
+ // ------------------------------------------------------------------------
+
+
+ // ----- Environment --------------------------------------------------
+ std::string srcDir = gSystem->Getenv("VMCWORKDIR"); // top source directory
+ // ------------------------------------------------------------------------
+
+ std::unique_ptr<CbmStsAnalysis> run = std::make_unique<CbmStsAnalysis>();
+
+ FairFileSource* inputSource = new FairFileSource(rec_file.c_str());
+ inputSource->AddFriend(raw_file.c_str());
+
+ FairRootFileSink* outputSink = new FairRootFileSink(out_file.c_str());
+
+
+ if (bAli) {
+ TString alignmentMatrixFileName = srcDir + "/parameters/mcbm/AlignmentMatrices_" + geo_setup_tag + ".root";
+ if (alignmentMatrixFileName.Length() != 0) {
+ std::map<std::string, TGeoHMatrix>* matrices{nullptr};
+ TFile* misalignmentMatrixRootfile = new TFile(alignmentMatrixFileName, "READ");
+ if (misalignmentMatrixRootfile->IsOpen()) {
+ gDirectory->GetObject("MisalignMatrices", matrices);
+ misalignmentMatrixRootfile->Close();
+ }
+ else {
+ LOG(error) << "Could not open file " << alignmentMatrixFileName << "\n Exiting";
+ exit(1);
+ }
+
+ if (matrices) {
+ run->AddAlignmentMatrices(*matrices);
+ }
+ else {
+ LOG(error) << "Alignment required but no matrices found."
+ << "\n Exiting";
+ exit(1);
+ }
+ }
+ }
+
+ run->SetGeomFile(geo_file.c_str());
+ run->SetSource(inputSource);
+ run->SetSink(outputSink);
+
+
+ // ------------------------------------------------------------------------
+ // Configure filters
+ // ------------------------------------------------------------------------
+ CbmCutMap ana_filter;
+ ana_filter.AddCbmCut(CbmCutId::kEventNofGlobalTrack)->SetMin(1);
+ ana_filter.AddCbmCut(CbmCutId::kGlobalTrackStsSize)->SetMin(2);
+ ana_filter.AddCbmCut(CbmCutId::kGlobalTrackTrdSize)->SetMin(1);
+ ana_filter.AddCbmCut(CbmCutId::kGlobalTrackTofSize)->SetMin(1);
+ ana_filter.AddCbmCut(CbmCutId::kGlobalTrackPval)->SetMax(0.5);
+
+ // ------------------------------------------------------------------------
+ // Configure analysis
+ // ------------------------------------------------------------------------
+ CbmStsResolution* sts_resolution = new CbmStsResolution(0.5);
+ sts_resolution->SetCutMap(&ana_filter);
+ sts_resolution->EnableRefSensor(0x10000002); // U0 L0 U0
+ sts_resolution->EnableRefSensor(0x10018422); // U1 L1 U1
+ sts_resolution->SetMinTrkHitQ(10000);
+ run->AddTask(sts_resolution);
+
+ run->Init();
+ run->Run(0, n_of_ts);
+
+ ana_filter.Print();
+ // ------------------------------------------------------------------------
+ timer.Stop();
+ Double_t rtime = timer.RealTime();
+ Double_t ctime = timer.CpuTime();
+ cout << endl << endl;
+ cout << "Macro finished successfully." << endl;
+ cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
+ cout << endl;
+}
diff --git a/macro/sts/sts_time.C b/macro/sts/sts_time.C
new file mode 100644
index 0000000000000000000000000000000000000000..64a4013a319c6d79eb64221caeda09e055b6a6af
--- /dev/null
+++ b/macro/sts/sts_time.C
@@ -0,0 +1,56 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dario Ramirez [committer] */
+
+void sts_time(int run_id = 2984, int n_of_ts = 2, std::string raw_file = "",
+ std::string out_file = "cbm_sts_time_cal.root")
+{
+ // ----- Timer --------------------------------------------------------
+ TStopwatch timer;
+ timer.Start();
+
+ // --- Logger settings ----------------------------------------------------
+ TString logLevel = "DEBUG";
+ TString logVerbosity = "veryhigh";
+ logLevel = "INFO";
+ logVerbosity = "low";
+
+ FairLogger::GetLogger()->SetLogScreenLevel(logLevel.Data());
+ FairLogger::GetLogger()->SetLogVerbosityLevel(logVerbosity.Data());
+ // ------------------------------------------------------------------------
+
+
+ // ----- Environment --------------------------------------------------
+ std::string srcDir = gSystem->Getenv("VMCWORKDIR"); // top source directory
+ // ------------------------------------------------------------------------
+
+ std::unique_ptr<CbmStsAnalysis> run = std::make_unique<CbmStsAnalysis>();
+ FairFileSource* inputSource = new FairFileSource(raw_file.c_str());
+ FairRootFileSink* outputSink = new FairRootFileSink(out_file.c_str());
+
+ run->SetSource(inputSource);
+ run->SetSink(outputSink);
+
+ // ------------------------------------------------------------------------
+ // Configure Analysis filters
+ // ------------------------------------------------------------------------
+ CbmCutMap ana_filter;
+ ana_filter.AddCbmCut(CbmCutId::kBmonDigiSide)->SetRange(1, 1);
+
+ CbmStsTimeCal* sts_time = new CbmStsTimeCal(run_id, ECbmModuleId::kBmon, -100, 100);
+ sts_time->SetCutMap(&ana_filter);
+ sts_time->SetWalkFile(Form("%s/parameters/online/StsWalkMap_mcbm2024.yaml", srcDir.c_str()));
+ run->AddTask(sts_time);
+
+ run->Init();
+ run->Run(0, n_of_ts);
+
+ // ------------------------------------------------------------------------
+ timer.Stop();
+ Double_t rtime = timer.RealTime();
+ Double_t ctime = timer.CpuTime();
+ cout << endl << endl;
+ cout << "Macro finished successfully." << endl;
+ cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
+ cout << endl;
+}