From b00009328c0e93436766d319ee64856981829c1a Mon Sep 17 00:00:00 2001
From: Evgeny Lavrik <evgeny.lavrik@uni-tuebingen.de>
Date: Mon, 7 Sep 2020 16:00:30 +0200
Subject: [PATCH] Add dEdx field to StsTrack, add dEdx calculation method to
CbmStsTrackFinder
---
core/data/sts/CbmStsTrack.h | 16 ++++-
reco/base/CMakeLists.txt | 4 ++
reco/base/CbmStsTrackFinder.cxx | 104 ++++++++++++++++++++++++++++++--
reco/base/CbmStsTrackFinder.h | 18 ++++++
4 files changed, 136 insertions(+), 6 deletions(-)
diff --git a/core/data/sts/CbmStsTrack.h b/core/data/sts/CbmStsTrack.h
index eadc7cf3f..737fb68ee 100644
--- a/core/data/sts/CbmStsTrack.h
+++ b/core/data/sts/CbmStsTrack.h
@@ -111,6 +111,18 @@ public:
virtual std::string ToString() const;
+ /** Get energy loss
+ ** @return median energy loss
+ **/
+ Float_t GetELoss() const { return fELoss; }
+
+ /** Set energy loss
+ ** @param median energy loss
+ **/
+ void SetELoss(Float_t ELoss) { fELoss = ELoss; }
+
+ constexpr static Float_t ELossOverflow() { return 1.e6; }
+
private:
/** Array with indices of the MVD hits attached to the track **/
std::vector<Int_t> fMvdHitIndex;
@@ -119,8 +131,10 @@ private:
/** Impact parameter of track at target z, in units of its error **/
Double32_t fB;
+ /** median dE/dx [e/300µm] **/
+ Float_t fELoss {-1.f};
- ClassDef(CbmStsTrack, 2);
+ ClassDef(CbmStsTrack, 3);
};
#endif
diff --git a/reco/base/CMakeLists.txt b/reco/base/CMakeLists.txt
index 7ec8faedf..43b1f1148 100644
--- a/reco/base/CMakeLists.txt
+++ b/reco/base/CMakeLists.txt
@@ -39,6 +39,10 @@ set(INCLUDE_DIRECTORIES
#${CBMDATA_DIR}
#${CBMROOT_SOURCE_DIR}/reco
${CBMROOT_SOURCE_DIR}/reco/base
+${CBMROOT_SOURCE_DIR}/core/base
+${CBMROOT_SOURCE_DIR}/core/data
+${CBMROOT_SOURCE_DIR}/core/data/base
+${CBMROOT_SOURCE_DIR}/core/data/sts
)
set(SYSTEM_INCLUDE_DIRECTORIES
diff --git a/reco/base/CbmStsTrackFinder.cxx b/reco/base/CbmStsTrackFinder.cxx
index 20783962a..dc64f6570 100644
--- a/reco/base/CbmStsTrackFinder.cxx
+++ b/reco/base/CbmStsTrackFinder.cxx
@@ -7,12 +7,14 @@
// Empty file, just there to please CINT
#include "CbmStsTrackFinder.h"
+#include "CbmDigiManager.h"
+#include "CbmStsCluster.h"
+#include "CbmStsDigi.h"
+#include "CbmStsHit.h"
+#include "CbmStsTrack.h"
+#include "FairRootManager.h"
-//#include "CbmStsDigiScheme.h"
-
-//#include "FairField.h"
-
-//#include "TClonesArray.h"
+#include "TClonesArray.h"
CbmStsTrackFinder::CbmStsTrackFinder()
: TNamed()
@@ -21,6 +23,98 @@ CbmStsTrackFinder::CbmStsTrackFinder()
, fMvdHits(nullptr)
, fStsHits(nullptr)
, fTracks(nullptr)
+ , fStsClusters(nullptr)
, fVerbose(0) {}
+double CbmStsTrackFinder::VecMedian(std::vector<double>& vec) {
+ if (vec.empty()) { return 0.; }
+
+ auto mid = vec.size() / 2;
+ std::nth_element(vec.begin(), vec.begin() + mid, vec.end());
+ auto median = vec[mid];
+ if (!(vec.size() & 1)) {
+ auto max_it = std::max_element(vec.begin(), vec.begin() + mid);
+ median = (*max_it + median) / 2.0;
+ }
+ return median;
+}
+
+double CbmStsTrackFinder::CalculateEloss(CbmStsTrack* cbmStsTrack) {
+ if (!fStsClusters) {
+ FairRootManager* ioman = FairRootManager::Instance();
+ assert(ioman);
+
+ fStsClusters = (TClonesArray*) ioman->GetObject("StsCluster");
+ assert(fStsClusters);
+ }
+
+ CbmDigiManager* digiManager = CbmDigiManager::Instance();
+
+ std::vector<double> dEdxAllveto;
+
+ double dr = 1.;
+ for (int iHit = 0; iHit < cbmStsTrack->GetNofStsHits(); ++iHit) {
+ bool frontVeto = kFALSE, backVeto = kFALSE;
+ CbmStsHit* stsHit =
+ (CbmStsHit*) fStsHits->At(cbmStsTrack->GetStsHitIndex(iHit));
+
+ double x, y, z, xNext, yNext, zNext;
+ x = stsHit->GetX();
+ y = stsHit->GetY();
+ z = stsHit->GetZ();
+
+ if (iHit != cbmStsTrack->GetNofStsHits() - 1) {
+ CbmStsHit* stsHitNext =
+ (CbmStsHit*) fStsHits->At(cbmStsTrack->GetStsHitIndex(iHit + 1));
+ xNext = stsHitNext->GetX();
+ yNext = stsHitNext->GetY();
+ zNext = stsHitNext->GetZ();
+ dr = sqrt((xNext - x) * (xNext - x) + (yNext - y) * (yNext - y)
+ + (zNext - z) * (zNext - z))
+ / (zNext - z); // if *300um, you get real reconstructed dr
+ } // else dr stay previous
+
+ CbmStsCluster* frontCluster =
+ (CbmStsCluster*) fStsClusters->At(stsHit->GetFrontClusterId());
+ CbmStsCluster* backCluster =
+ (CbmStsCluster*) fStsClusters->At(stsHit->GetBackClusterId());
+
+ if (!frontCluster || !backCluster) {
+ LOG(info)
+ << "CbmStsTrackFinder::CalculateEloss: no front or back cluster";
+ continue;
+ }
+
+ //check if at least one digi in a cluster has overflow --- charge is registered in the last ADC channel #31
+ for (int iDigi = 0; iDigi < frontCluster->GetNofDigis(); ++iDigi) {
+ if (CbmStsTrackFinder::MaxAdcVal()
+ == (digiManager->Get<CbmStsDigi>(frontCluster->GetDigi(iDigi))
+ ->GetCharge()))
+ frontVeto = kTRUE;
+ }
+ for (int iDigi = 0; iDigi < backCluster->GetNofDigis(); ++iDigi) {
+ if (CbmStsTrackFinder::MaxAdcVal()
+ == (digiManager->Get<CbmStsDigi>(backCluster->GetDigi(iDigi))
+ ->GetCharge()))
+ backVeto = kTRUE;
+ }
+
+ if (!frontVeto) dEdxAllveto.push_back((frontCluster->GetCharge()) / dr);
+ if (!backVeto) dEdxAllveto.push_back((backCluster->GetCharge()) / dr);
+ }
+
+ float dEdXSTS = CbmStsTrack::ELossOverflow();
+ if (dEdxAllveto.size() != 0) dEdXSTS = VecMedian(dEdxAllveto);
+ return dEdXSTS;
+}
+
+void CbmStsTrackFinder::FillEloss() {
+ Int_t nStsTracks = fTracks->GetEntriesFast();
+ for (Int_t stsTrackIndex = 0; stsTrackIndex < nStsTracks; stsTrackIndex++) {
+ CbmStsTrack* cbmStsTrack = (CbmStsTrack*) fTracks->At(stsTrackIndex);
+ double dEdXSTS = CalculateEloss(cbmStsTrack);
+ cbmStsTrack->SetELoss(dEdXSTS);
+ }
+}
+
ClassImp(CbmStsTrackFinder)
diff --git a/reco/base/CbmStsTrackFinder.h b/reco/base/CbmStsTrackFinder.h
index 34ced83bf..72ab76d68 100644
--- a/reco/base/CbmStsTrackFinder.h
+++ b/reco/base/CbmStsTrackFinder.h
@@ -24,6 +24,7 @@ class TClonesArray;
class CbmStsDigiScheme;
class FairField;
class CbmEvent;
+class CbmStsTrack;
class CbmStsTrackFinder : public TNamed {
@@ -80,9 +81,26 @@ protected:
TClonesArray* fMvdHits; // MvdHit array
TClonesArray* fStsHits; // StsHit array
TClonesArray* fTracks; // StsTrack array
+ TClonesArray* fStsClusters; // StsCluster array
Int_t fVerbose; // Verbosity level
+
+ /** Median energy loss calculation for the tracks in event/timeslice
+ ** Ported from CbmKFParticleFinderPID
+ ** Description of the method given at 30th CBM CM
+ ** https://indico.gsi.de/event/4760/session/4/contribution/80/material/slides/0.pdf
+ **/
+ double CalculateEloss(CbmStsTrack* cbmStsTrack);
+ void FillEloss();
+
private:
+ constexpr static int MaxAdcVal() { return 31; }
+
+ /** Calculate median value of a vector
+ **/
+ double VecMedian(std::vector<double>& vec);
+
+
CbmStsTrackFinder(const CbmStsTrackFinder&);
CbmStsTrackFinder& operator=(const CbmStsTrackFinder&);
--
GitLab