diff --git a/core/data/sts/CbmStsTrack.h b/core/data/sts/CbmStsTrack.h index eadc7cf3fdef003deb873c26661ebdaf815377d1..737fb68eec8b2f074a943427845d4ea6f347e6ac 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 7ec8faedf81e7f19cd8c8255c0f14c8edf4cd5cf..43b1f114849deea53a6d7385a68d474dc60cc72e 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 20783962a66738d7c5300aa71665d0a2980ca339..dc64f6570b3846e48bac6f6ec0856dd3260023cc 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 34ced83bf81312d6dfc3306e5dad5942a5d5ef6b..72ab76d68abebe8754eae97450d3a9285d7a9a11 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&);