Add dEdx field to StsTrack, add dEdx calculation method to CbmStsTrackFinder

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

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

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)

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&);