From 3d7023d33643fc91d6e07d44b82e48212ff1bd50 Mon Sep 17 00:00:00 2001 From: RealAurio Date: Mon, 21 Jul 2025 17:48:44 +0200 Subject: [PATCH 1/9] Refactored a bunch of the code related to the fitting step of the analysis. --- AtData/AtDataLinkDef.h | 7 + AtData/AtFitResult.cxx | 3 + AtData/AtFitResult.h | 49 ++++ AtData/AtFitTrackResult.cxx | 16 ++ AtData/AtFitTrackResult.h | 46 +++ AtData/AtFittedTrack.cxx | 60 ++-- AtData/AtFittedTrack.h | 183 ++++++------ AtData/AtFittedTrackOld.cxx | 45 +++ AtData/AtFittedTrackOld.h | 144 ++++++++++ AtData/AtMCResult.cxx | 8 +- AtData/AtMCResult.h | 10 +- AtData/AtMCResultOld.cxx | 14 + AtData/AtMCResultOld.h | 35 +++ AtData/AtTrackingEventOld.cxx | 48 ++++ AtData/AtTrackingEventOld.h | 59 ++++ AtData/CMakeLists.txt | 12 +- AtReconstruction/AtFitter/AtFitter.cxx | 267 +---------------- AtReconstruction/AtFitter/AtFitter.h | 62 ++-- AtReconstruction/AtFitter/AtFitterOld.cxx | 269 ++++++++++++++++++ AtReconstruction/AtFitter/AtFitterOld.h | 53 ++++ AtReconstruction/AtFitter/AtMCFitter.cxx | 4 +- AtReconstruction/AtFitter/AtMCFitterOld.cxx | 216 ++++++++++++++ AtReconstruction/AtFitter/AtMCFitterOld.h | 136 +++++++++ .../AtFitter/AtMCFitterTaskOld.cxx | 51 ++++ AtReconstruction/AtFitter/AtMCFitterTaskOld.h | 43 +++ AtReconstruction/AtFitterTask.cxx | 65 ++++- AtReconstruction/AtFitterTask.h | 27 +- AtReconstruction/AtFitterTaskOld.cxx | 108 +++++++ AtReconstruction/AtFitterTaskOld.h | 79 +++++ AtReconstruction/AtReconstructionLinkDef.h | 11 +- AtReconstruction/CMakeLists.txt | 13 +- AtTools/AtKinematics.cxx | 9 + AtTools/AtKinematics.h | 1 + AtTools/AtTrackTransformer.cxx | 132 +++++++++ AtTools/AtTrackTransformer.h | 5 + 35 files changed, 1854 insertions(+), 436 deletions(-) create mode 100644 AtData/AtFitResult.cxx create mode 100644 AtData/AtFitResult.h create mode 100644 AtData/AtFitTrackResult.cxx create mode 100644 AtData/AtFitTrackResult.h create mode 100644 AtData/AtFittedTrackOld.cxx create mode 100644 AtData/AtFittedTrackOld.h create mode 100644 AtData/AtMCResultOld.cxx create mode 100644 AtData/AtMCResultOld.h create mode 100644 AtData/AtTrackingEventOld.cxx create mode 100644 AtData/AtTrackingEventOld.h create mode 100644 AtReconstruction/AtFitter/AtFitterOld.cxx create mode 100644 AtReconstruction/AtFitter/AtFitterOld.h create mode 100644 AtReconstruction/AtFitter/AtMCFitterOld.cxx create mode 100644 AtReconstruction/AtFitter/AtMCFitterOld.h create mode 100644 AtReconstruction/AtFitter/AtMCFitterTaskOld.cxx create mode 100644 AtReconstruction/AtFitter/AtMCFitterTaskOld.h create mode 100644 AtReconstruction/AtFitterTaskOld.cxx create mode 100644 AtReconstruction/AtFitterTaskOld.h diff --git a/AtData/AtDataLinkDef.h b/AtData/AtDataLinkDef.h index b1cc7f9da..013860b29 100644 --- a/AtData/AtDataLinkDef.h +++ b/AtData/AtDataLinkDef.h @@ -43,5 +43,12 @@ #pragma link C++ enum AtPatterns::PatternType; #pragma link C++ function AtPatterns::CreatePattern; +#pragma link C++ class AtFitResult + ; +#pragma link C++ class AtFitTrackResult + ; #pragma link C++ class MCFitter::AtMCResult + ; + +#pragma link C++ class AtTrackingEventOld + ; +#pragma link C++ class AtFittedTrackOld + ; +#pragma link C++ class MCFitter::AtMCResultOld + ; + #endif diff --git a/AtData/AtFitResult.cxx b/AtData/AtFitResult.cxx new file mode 100644 index 000000000..5412b7864 --- /dev/null +++ b/AtData/AtFitResult.cxx @@ -0,0 +1,3 @@ +#include "AtFitResult.h" + +ClassImp(AtFitResult); diff --git a/AtData/AtFitResult.h b/AtData/AtFitResult.h new file mode 100644 index 000000000..5542d29ff --- /dev/null +++ b/AtData/AtFitResult.h @@ -0,0 +1,49 @@ +#ifndef ATFITRESULT_H +#define ATFITRESULT_H + +#include "AtFitTrackResult.h" + +#include + +#include // for Double_t, THashConsistencyHolder, ClassDefOverride +#include + +#include +#include +#include +#include + +class TBuffer; +class TClass; +class TMemberInspector; + +/** + * Class for storing the result of the fit for the entire AtTrackingEvent from an AtFitter class. + */ +class AtFitResult : public TObject { +public: + using TrackResultPtr = std::unique_ptr; + using TrackResultsVector = std::vector; + using ResultsMap = std::map; + +protected: + // Vector to store the results for all different fits done to all tracks in the event. + ResultsMap fResults; + + // Event ID for which this fit was done. + ULong_t fEventID; + +public: + AtFitResult() = default; + ~AtFitResult() = default; + + void SetTrackResultsVector(Int_t trackID, TrackResultsVector results) { fResults[trackID] = std::move(results); } + + void SetEventID(ULong_t id) { fEventID = id; } + + TrackResultsVector &GetTrackResultsVector(Int_t trackID) { return fResults.at(trackID); } + + ClassDefOverride(AtFitResult, 1); +}; + +#endif diff --git a/AtData/AtFitTrackResult.cxx b/AtData/AtFitTrackResult.cxx new file mode 100644 index 000000000..2c93509b9 --- /dev/null +++ b/AtData/AtFitTrackResult.cxx @@ -0,0 +1,16 @@ +#include "AtFitTrackResult.h" + +#include +ClassImp(AtFitTrackResult); + +void AtFitTrackResult::Print() const + +{ + std::cout << " Fit result for track with ID " << fTrackID << ":" << std::endl; + + std::cout << " Statistics: " << std::endl; + std::cout << " PValue = " << fPValue << std::endl; + std::cout << " Chi2 = " << fChi2 << std::endl; + std::cout << " NDF = " << fNdf << std::endl; + std::cout << " Converged = " << fFitConverged << std::endl; +} diff --git a/AtData/AtFitTrackResult.h b/AtData/AtFitTrackResult.h new file mode 100644 index 000000000..2960a84c9 --- /dev/null +++ b/AtData/AtFitTrackResult.h @@ -0,0 +1,46 @@ +#ifndef ATFITTRACKRESULT_H +#define ATFITTRACKRESULT_H + +#include // for Double_t, THashConsistencyHolder, ClassDefOverride +#include + +class TBuffer; +class TClass; +class TMemberInspector; + +/** + * Class for storing the result of the fit of an AtTrack from an AtFitter class. + */ +class AtFitTrackResult : public TObject { +protected: + // Copy of the statistics parameters which are also saved in AtFittedTrack. + Double_t fPValue{0}; + Double_t fChi2{0}; + Int_t fNdf{0}; + Bool_t fFitConverged{false}; + + // The track ID for which this fit was done for. + Int_t fTrackID{-1}; + +public: + AtFitTrackResult() = default; + ~AtFitTrackResult() = default; + + void SetPValue(Double_t value) { fPValue = value; } + void SetChi2(Double_t value) { fChi2 = value; } + void SetNdf(Int_t value) { fNdf = value; } + void SetFitConverged(Bool_t value) { fFitConverged = value; } + void SetTrackID(Int_t value) { fTrackID = value; } + + Double_t GetPValue() const { return fPValue; } + Double_t GetChi2() const { return fChi2; } + Int_t GetNdf() const { return fNdf; } + Bool_t GetFitConverged() const { return fFitConverged; } + Int_t GetTrackID() const { return fTrackID; } + + virtual void Print() const; + + ClassDefOverride(AtFitTrackResult, 1); +}; + +#endif diff --git a/AtData/AtFittedTrack.cxx b/AtData/AtFittedTrack.cxx index 4481bf8f7..9a5756bc2 100644 --- a/AtData/AtFittedTrack.cxx +++ b/AtData/AtFittedTrack.cxx @@ -1,45 +1,41 @@ #include "AtFittedTrack.h" -#include - -#include -#include - ClassImp(AtFittedTrack); -using XYZVector = ROOT::Math::XYZVector; - -const std::tuple AtFittedTrack::GetEnergyAngles() +void AtFittedTrack::SetKinematics(int particleIdx, Double_t energy, Double_t theta, Double_t phi, Double_t energyxtr, + Double_t thetaxtr, Double_t phixtr) { - return std::forward_as_tuple(fEnergy, fEnergyXtr, fTheta, fPhi, fEnergyPRA, fThetaPRA, fPhiPRA); + while (particleIdx >= fKinematics.size()) { + Kinematics newKinematics; + fKinematics.push_back(newKinematics); + } + + fKinematics[particleIdx].kineticEnergy = energy; + fKinematics[particleIdx].theta = theta; + fKinematics[particleIdx].phi = phi; + fKinematics[particleIdx].kineticEnergyXtr = energyxtr; + fKinematics[particleIdx].thetaXtr = thetaxtr; + fKinematics[particleIdx].phiXtr = phixtr; } -const std::tuple AtFittedTrack::GetVertices() +void AtFittedTrack::SetParticleInfo(int particleIdx, std::string pdg, Int_t charge, Double_t mass) { - return std::forward_as_tuple(fInitialPos, fInitialPosPRA, fInitialPosXtr); + while (particleIdx >= fParticleInfo.size()) { + ParticleInfo newParticleInfo; + fParticleInfo.push_back(newParticleInfo); + } + + fParticleInfo[particleIdx].idPDG = TString(pdg); + fParticleInfo[particleIdx].charge = charge; + fParticleInfo[particleIdx].mass = mass; } -const std::tuple AtFittedTrack::GetStats() +void AtFittedTrack::SetVertex(int particleIdx, XYZVector point) { - return std::forward_as_tuple(fPValue, fChi2, fBChi2, fNdf, fBNdf, fFitConverged); -} + while (particleIdx >= fVertex.size()) { + XYZVector newVertex; + fVertex.push_back(newVertex); + } -const std::tuple AtFittedTrack::GetTrackProperties() -{ - return std::forward_as_tuple(fCharge, fBrho, fELossADC, fDEdxADC, fPDG, fTrackPoints); + fVertex[particleIdx] = point; } - -const std::tuple AtFittedTrack::GetIonChamber() -{ - return std::forward_as_tuple(fIonChamberEnergy, fIonChamberTime); -} - -const std::tuple AtFittedTrack::GetExcitationEnergy() -{ - return std::forward_as_tuple(fExcitationEnergy, fExcitationEnergyXtr); -} - -const std::tuple AtFittedTrack::GetDistances() -{ - return std::forward_as_tuple(fDistanceXtr, fTrackLength, fPOCAXtr); -} \ No newline at end of file diff --git a/AtData/AtFittedTrack.h b/AtData/AtFittedTrack.h index 4cb0c1044..fce53b138 100644 --- a/AtData/AtFittedTrack.h +++ b/AtData/AtFittedTrack.h @@ -1,6 +1,8 @@ #ifndef ATFITTEDTRACK_H #define ATFITTEDTRACK_H +#include "AtFitTrackResult.h" + #include #include #include @@ -8,6 +10,7 @@ #include #include #include +#include #include #include @@ -20,125 +23,129 @@ class TClass; class TMemberInspector; class AtFittedTrack : public TObject { - +public: using XYZVector = ROOT::Math::XYZVector; + using TrackResultPtr = std::unique_ptr; + + struct Kinematics { + Double_t kineticEnergy{-1}; // Kinetic energy + Double_t theta{-1}; // Theta scattering angle + Double_t phi{-1}; // Phi scattering angle + Double_t kineticEnergyXtr{-1}; // Extrapolated kinetic energy + Double_t thetaXtr{-1}; // Extrapolated theta scattering angle + Double_t phiXtr{-1}; // Extrapolated phi scattering angle + }; + + struct Statistics { + Double_t pValue{-1}; // Probability for rejecting the fit hypothesis + Double_t chi2{-1}; // Chi2 of the fit + Int_t NDF{-1}; // Number of degrees of freedom for the fit + Double_t redChi2{-1}; // Reduced chi2 + Bool_t fitConverged{0}; // Whether or not the fit managed to converge + }; + + struct ParticleInfo { + TString idPDG{""}; // PDG code of the particle + Int_t charge{0}; // Charge of the particle + Double_t mass{-1}; // Mass of the particle + }; + + struct TrackProperties { + XYZVector initialPosition; // Position of the first hit + XYZVector initialPositionXtr; // Position of the point closest to (0,0) + Double_t extrapolatedDistance{-1}; // Distance initialPosition->initialPositionXtr along the pattern + Double_t distancePOCA{-1}; // Distance initialPositionXtr->(0,0) + Double_t trackLength{-1}; // Distance initialPosition->End of charge + Double_t trackLengthXtr{-1}; // Distance initialPositionXtr->End of charge + Double_t estimateTotalCharge{-1}; // Sum of the charge of all hits + Double_t estimateDeDx{-1}; // Sum of the charge of all hits divided by range + Int_t trackPoints{-1}; // Number of hits in the track + }; private: Int_t fTrackID{-1}; //< Track ID from pattern recognition - Float_t fEnergy{0}; - Float_t fTheta{0}; - Float_t fPhi{0}; - Float_t fEnergyPRA{0}; - Float_t fThetaPRA{0}; - Float_t fPhiPRA{0}; - - Float_t fExcitationEnergy{0}; + // Kinematic variables obtained by the fit. + std::vector fKinematics; - XYZVector fInitialPos; // xiniFitVec,yiniFitVec,ziniFitVec; - XYZVector fInitialPosPRA; // xiniPRAVec,yiniPRAVec,ziniPRAVec; - XYZVector fInitialPosXtr; + // Particle information. + std::vector fParticleInfo; - Float_t fIonChamberEnergy{0}; - Int_t fIonChamberTime{0}; + // Vertex where the particle has originated from. + std::vector fVertex; - Float_t fEnergyXtr{0}; - Float_t fExcitationEnergyXtr{0}; + // Track properties. + TrackProperties fTrackProperties; - Float_t fDistanceXtr{0}; - Float_t fTrackLength{0}; - Float_t fPOCAXtr{0}; + // Parameters regarding the statistics of the fit. + Statistics fStats; - Float_t fPValue{0}; - Float_t fChi2{0}; - Float_t fBChi2{0}; - Float_t fNdf{0}; - Float_t fBNdf{0}; - Bool_t fFitConverged{0}; - - Int_t fCharge{0}; - Float_t fBrho{0}; - Float_t fELossADC{0}; - Float_t fDEdxADC{0}; - std::string fPDG{0}; - Int_t fTrackPoints{0}; + // Copy of the AtFitTrackResult object corresponding to the fit used for this track. + TrackResultPtr fTrackResult{nullptr}; public: AtFittedTrack() = default; ~AtFittedTrack() = default; - inline void SetTrackID(Int_t trackid) { fTrackID = trackid; } + void SetTrackID(Int_t trackid) { fTrackID = trackid; } - inline void SetEnergyAngles(Float_t energy, Float_t energyxtr, Float_t theta, Float_t phi, Float_t energypra, - Float_t thetapra, Float_t phipra) - { - fEnergy = energy; - fEnergyXtr = energyxtr; - fTheta = theta; - fPhi = phi; - fEnergyPRA = energypra; - fThetaPRA = thetapra; - fPhiPRA = phi; - } + void SetKinematics(int particleIdx, Double_t energy, Double_t theta, Double_t phi, Double_t energyxtr, + Double_t thetaxtr, Double_t phixtr); + void SetParticleInfo(int particleIdx, std::string pdg, Int_t charge, Double_t mass); + void SetVertex(int particleIdx, XYZVector point); - inline void SetVertexPosition(XYZVector inipos, XYZVector iniposPRA, XYZVector iniposxtr) + void + SetKinematics(Double_t energy, Double_t theta, Double_t phi, Double_t energyxtr, Double_t thetaxtr, Double_t phixtr) { - fInitialPos = inipos; - fInitialPosPRA = iniposPRA; - fInitialPosXtr = iniposxtr; + SetKinematics(0, energy, theta, phi, energyxtr, thetaxtr, phixtr); } - inline void SetStats(Float_t pvalue, Float_t chi2, Float_t bchi2, Float_t ndf, Float_t bndf, Bool_t conv) - { - fPValue = pvalue; - fChi2 = chi2; - fBChi2 = bchi2; - fNdf = ndf; - fBNdf = bndf; - fFitConverged = conv; - } + void SetParticleInfo(std::string pdg, Int_t charge, Double_t mass) { SetParticleInfo(0, pdg, charge, mass); } - inline void - SetTrackProperties(Int_t charge, Float_t brho, Float_t eloss, Float_t dedx, std::string pdg, Int_t points) - { - fCharge = charge; - fBrho = brho; - fELossADC = eloss; - fDEdxADC = dedx; - fPDG = pdg; - fTrackPoints = points; - } + void SetVertex(XYZVector point) { SetVertex(0, point); } - inline void SetIonChamber(Float_t icenergy, Int_t ictime) + void SetTrackProperties(XYZVector initialPosition, XYZVector initialPositionXtr, Double_t extrapolatedDistance, + Double_t distancePOCA, Double_t trackLength, Double_t trackLengthXtr, + Double_t estimateTotalCharge, Int_t trackPoints) { - fIonChamberEnergy = icenergy; - fIonChamberTime = ictime; + fTrackProperties.initialPosition = initialPosition; + fTrackProperties.initialPositionXtr = initialPositionXtr; + fTrackProperties.extrapolatedDistance = extrapolatedDistance; + fTrackProperties.distancePOCA = distancePOCA; + fTrackProperties.trackLength = trackLength; + fTrackProperties.trackLengthXtr = trackLengthXtr; + fTrackProperties.estimateTotalCharge = estimateTotalCharge; + fTrackProperties.estimateDeDx = estimateTotalCharge / trackLengthXtr; + fTrackProperties.trackPoints = trackPoints; } - inline void SetExcitationEnergy(Float_t exenergy, Float_t exenergyxtr) + void SetStats(Double_t pvalue, Double_t chi2, Int_t ndf, Bool_t conv) { - fExcitationEnergy = exenergy; - fExcitationEnergyXtr = exenergyxtr; + fStats.pValue = pvalue; + fStats.chi2 = chi2; + fStats.NDF = ndf; + fStats.redChi2 = chi2 / ndf; + fStats.fitConverged = conv; } - inline void SetDistances(Float_t distancextr, Float_t length, Float_t poca) - { - fDistanceXtr = distancextr; - fTrackLength = length; - fPOCAXtr = poca; - } + void SetTrackResult(TrackResultPtr trackResult) { fTrackResult = std::move(trackResult); } const Int_t GetTrackID() { return fTrackID; } - const std::tuple GetEnergyAngles(); - const std::tuple GetVertices(); - const std::tuple GetStats(); - const std::tuple GetTrackProperties(); - const std::tuple GetIonChamber(); - const std::tuple GetExcitationEnergy(); - const std::tuple GetDistances(); + const Kinematics GetKinematics(int particleIdx) { return fKinematics[particleIdx]; } + const ParticleInfo GetParticleInfo(int particleIdx) { return fParticleInfo[particleIdx]; } + const XYZVector GetVertex(int particleIdx) { return fVertex[particleIdx]; } + + const Kinematics GetKinematics() { return fKinematics[0]; } + const ParticleInfo GetParticleInfo() { return fParticleInfo[0]; } + const XYZVector GetVertex() { return fVertex[0]; } + + const TrackProperties GetTrackProperties() { return fTrackProperties; } + const Statistics GetStats() { return fStats; } + + TrackResultPtr &GetTrackResult() { return fTrackResult; } - ClassDef(AtFittedTrack, 1); + ClassDef(AtFittedTrack, 2); }; -#endif \ No newline at end of file +#endif diff --git a/AtData/AtFittedTrackOld.cxx b/AtData/AtFittedTrackOld.cxx new file mode 100644 index 000000000..8b222a55e --- /dev/null +++ b/AtData/AtFittedTrackOld.cxx @@ -0,0 +1,45 @@ +#include "AtFittedTrackOld.h" + +#include + +#include +#include + +ClassImp(AtFittedTrackOld); + +using XYZVector = ROOT::Math::XYZVector; + +const std::tuple AtFittedTrackOld::GetEnergyAngles() +{ + return std::forward_as_tuple(fEnergy, fEnergyXtr, fTheta, fPhi, fEnergyPRA, fThetaPRA, fPhiPRA); +} + +const std::tuple AtFittedTrackOld::GetVertices() +{ + return std::forward_as_tuple(fInitialPos, fInitialPosPRA, fInitialPosXtr); +} + +const std::tuple AtFittedTrackOld::GetStats() +{ + return std::forward_as_tuple(fPValue, fChi2, fBChi2, fNdf, fBNdf, fFitConverged); +} + +const std::tuple AtFittedTrackOld::GetTrackProperties() +{ + return std::forward_as_tuple(fCharge, fBrho, fELossADC, fDEdxADC, fPDG, fTrackPoints); +} + +const std::tuple AtFittedTrackOld::GetIonChamber() +{ + return std::forward_as_tuple(fIonChamberEnergy, fIonChamberTime); +} + +const std::tuple AtFittedTrackOld::GetExcitationEnergy() +{ + return std::forward_as_tuple(fExcitationEnergy, fExcitationEnergyXtr); +} + +const std::tuple AtFittedTrackOld::GetDistances() +{ + return std::forward_as_tuple(fDistanceXtr, fTrackLength, fPOCAXtr); +} diff --git a/AtData/AtFittedTrackOld.h b/AtData/AtFittedTrackOld.h new file mode 100644 index 000000000..4ef1ffcfb --- /dev/null +++ b/AtData/AtFittedTrackOld.h @@ -0,0 +1,144 @@ +#ifndef ATFITTEDTRACKOLD_H +#define ATFITTEDTRACKOLD_H + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +class TBuffer; +class TClass; +class TMemberInspector; + +class [[deprecated]] AtFittedTrackOld : public TObject { + + using XYZVector = ROOT::Math::XYZVector; + +private: + Int_t fTrackID{-1}; //< Track ID from pattern recognition + + Float_t fEnergy{0}; + Float_t fTheta{0}; + Float_t fPhi{0}; + Float_t fEnergyPRA{0}; + Float_t fThetaPRA{0}; + Float_t fPhiPRA{0}; + + Float_t fExcitationEnergy{0}; + + XYZVector fInitialPos; // xiniFitVec,yiniFitVec,ziniFitVec; + XYZVector fInitialPosPRA; // xiniPRAVec,yiniPRAVec,ziniPRAVec; + XYZVector fInitialPosXtr; + + Float_t fIonChamberEnergy{0}; + Int_t fIonChamberTime{0}; + + Float_t fEnergyXtr{0}; + Float_t fExcitationEnergyXtr{0}; + + Float_t fDistanceXtr{0}; + Float_t fTrackLength{0}; + Float_t fPOCAXtr{0}; + + Float_t fPValue{0}; + Float_t fChi2{0}; + Float_t fBChi2{0}; + Float_t fNdf{0}; + Float_t fBNdf{0}; + Bool_t fFitConverged{0}; + + Int_t fCharge{0}; + Float_t fBrho{0}; + Float_t fELossADC{0}; + Float_t fDEdxADC{0}; + std::string fPDG{0}; + Int_t fTrackPoints{0}; + +public: + AtFittedTrackOld() = default; + ~AtFittedTrackOld() = default; + + inline void SetTrackID(Int_t trackid) { fTrackID = trackid; } + + inline void SetEnergyAngles(Float_t energy, Float_t energyxtr, Float_t theta, Float_t phi, Float_t energypra, + Float_t thetapra, Float_t phipra) + { + fEnergy = energy; + fEnergyXtr = energyxtr; + fTheta = theta; + fPhi = phi; + fEnergyPRA = energypra; + fThetaPRA = thetapra; + fPhiPRA = phi; + } + + inline void SetVertexPosition(XYZVector inipos, XYZVector iniposPRA, XYZVector iniposxtr) + { + fInitialPos = inipos; + fInitialPosPRA = iniposPRA; + fInitialPosXtr = iniposxtr; + } + + inline void SetStats(Float_t pvalue, Float_t chi2, Float_t bchi2, Float_t ndf, Float_t bndf, Bool_t conv) + { + fPValue = pvalue; + fChi2 = chi2; + fBChi2 = bchi2; + fNdf = ndf; + fBNdf = bndf; + fFitConverged = conv; + } + + inline void + SetTrackProperties(Int_t charge, Float_t brho, Float_t eloss, Float_t dedx, std::string pdg, Int_t points) + { + fCharge = charge; + fBrho = brho; + fELossADC = eloss; + fDEdxADC = dedx; + fPDG = pdg; + fTrackPoints = points; + } + + inline void SetIonChamber(Float_t icenergy, Int_t ictime) + { + fIonChamberEnergy = icenergy; + fIonChamberTime = ictime; + } + + inline void SetExcitationEnergy(Float_t exenergy, Float_t exenergyxtr) + { + fExcitationEnergy = exenergy; + fExcitationEnergyXtr = exenergyxtr; + } + + inline void SetDistances(Float_t distancextr, Float_t length, Float_t poca) + { + fDistanceXtr = distancextr; + fTrackLength = length; + fPOCAXtr = poca; + } + + const Int_t GetTrackID() { return fTrackID; } + + const std::tuple GetEnergyAngles(); + const std::tuple GetVertices(); + const std::tuple GetStats(); + const std::tuple GetTrackProperties(); + const std::tuple GetIonChamber(); + const std::tuple GetExcitationEnergy(); + const std::tuple GetDistances(); + + ClassDef(AtFittedTrackOld, 1); +}; + +#endif diff --git a/AtData/AtMCResult.cxx b/AtData/AtMCResult.cxx index 23884150b..053e926fc 100644 --- a/AtData/AtMCResult.cxx +++ b/AtData/AtMCResult.cxx @@ -7,8 +7,12 @@ namespace MCFitter { void AtMCResult::Print() const { - std::cout << "Objective: " << fObjective << " Iteration: " << fIterNum << std::endl; + AtFitTrackResult::Print(); + std::cout << " MC fit specifics: " << std::endl; + + std::cout << " Iteration = " << fIterNum << std::endl; + std::cout << " Parameters:" << std::endl; for (auto &[name, val] : fParameters) - std::cout << name << ": " << val << std::endl; + std::cout << " " << name << " = " << val << std::endl; } } // namespace MCFitter diff --git a/AtData/AtMCResult.h b/AtData/AtMCResult.h index 74ea5d9f3..d4d48246b 100644 --- a/AtData/AtMCResult.h +++ b/AtData/AtMCResult.h @@ -1,6 +1,8 @@ #ifndef ATMCRESULT_H #define ATMCRESULT_H +#include "AtFitTrackResult.h" + #include // for Double_t, THashConsistencyHolder, ClassDefOverride #include @@ -15,19 +17,19 @@ namespace MCFitter { /** * Class for storing the result of an iteration in the AtMCFitter method. */ -class AtMCResult : public TObject { +class AtMCResult : public AtFitTrackResult { public: using ParamMap = std::map; - Double_t fObjective; //< Value f the objective function for this iteration ParamMap fParameters; //< Parameters used in simulation Int_t fIterNum; //< Iteration number. Used to map with the simulated event ID in the TTree. AtMCResult() = default; + ~AtMCResult() = default; - void Print() const; + void Print() const override; - ClassDefOverride(AtMCResult, 1); + ClassDefOverride(AtMCResult, 2); }; } // namespace MCFitter diff --git a/AtData/AtMCResultOld.cxx b/AtData/AtMCResultOld.cxx new file mode 100644 index 000000000..54481b7d4 --- /dev/null +++ b/AtData/AtMCResultOld.cxx @@ -0,0 +1,14 @@ +#include "AtMCResultOld.h" + +#include +ClassImp(MCFitter::AtMCResultOld); +namespace MCFitter { + +void AtMCResultOld::Print() const + +{ + std::cout << "Objective: " << fObjective << " Iteration: " << fIterNum << std::endl; + for (auto &[name, val] : fParameters) + std::cout << name << ": " << val << std::endl; +} +} // namespace MCFitter diff --git a/AtData/AtMCResultOld.h b/AtData/AtMCResultOld.h new file mode 100644 index 000000000..e095c7497 --- /dev/null +++ b/AtData/AtMCResultOld.h @@ -0,0 +1,35 @@ +#ifndef ATMCRESULTOLD_H +#define ATMCRESULTOLD_H + +#include // for Double_t, THashConsistencyHolder, ClassDefOverride +#include + +#include +#include // for string +class TBuffer; +class TClass; +class TMemberInspector; + +namespace MCFitter { + +/** + * Class for storing the result of an iteration in the AtMCFitter method. + */ +class [[deprecated]] AtMCResultOld : public TObject { +public: + using ParamMap = std::map; + + Double_t fObjective; //< Value f the objective function for this iteration + ParamMap fParameters; //< Parameters used in simulation + Int_t fIterNum; //< Iteration number. Used to map with the simulated event ID in the TTree. + + AtMCResultOld() = default; + + void Print() const; + + ClassDefOverride(AtMCResultOld, 1); +}; + +} // namespace MCFitter + +#endif // #ifndef ATMCRESULT_H diff --git a/AtData/AtTrackingEventOld.cxx b/AtData/AtTrackingEventOld.cxx new file mode 100644 index 000000000..8da4cc097 --- /dev/null +++ b/AtData/AtTrackingEventOld.cxx @@ -0,0 +1,48 @@ +#include "AtTrackingEventOld.h" + +#include +#include + +#include + +ClassImp(AtTrackingEventOld); + +AtTrackingEventOld::AtTrackingEventOld() : AtBaseEvent("Tracking Event") {} + +void AtTrackingEventOld::SetTrackArray(std::vector *trackArray) +{ + fTrackArray = *trackArray; +} +void AtTrackingEventOld::SetTrack(AtTrack *track) +{ + fTrackArray.push_back(*track); +} +void AtTrackingEventOld::SetVertex(Double_t vertex) +{ + fVertex = vertex; +} +void AtTrackingEventOld::SetGeoVertex(TVector3 vertex) +{ + fGeoVertex = vertex; +} +void AtTrackingEventOld::SetVertexEnergy(Double_t vertexEner) +{ + fVertexEnergy = vertexEner; +} + +/*std::vector AtTrackingEventOld::GetTrackArray() +{ + return fTrackArray; +}*/ +Double_t AtTrackingEventOld::GetVertex() +{ + return fVertex; +} +Double_t AtTrackingEventOld::GetVertexEnergy() +{ + return fVertexEnergy; +} +TVector3 AtTrackingEventOld::GetGeoVertex() +{ + return fGeoVertex; +} diff --git a/AtData/AtTrackingEventOld.h b/AtData/AtTrackingEventOld.h new file mode 100644 index 000000000..599f71b44 --- /dev/null +++ b/AtData/AtTrackingEventOld.h @@ -0,0 +1,59 @@ +#ifndef AtTRACKINGEVENTOLD_H +#define AtTRACKINGEVENTOLD_H + +#include "AtBaseEvent.h" +#include "AtFittedTrackOld.h" +#include "AtTrack.h" + +#include +#include +#include + +#include + +class TBuffer; +class TClass; +class TMemberInspector; + +class [[deprecated]] AtTrackingEventOld : public AtBaseEvent { + + using FTrackPtr = std::unique_ptr; + using FTrackVector = std::vector; + +public: + AtTrackingEventOld(); + virtual ~AtTrackingEventOld() = default; + + void SetTrackArray(std::vector *trackArray); + void SetTrack(AtTrack *track); + void SetVertex(Double_t vertex); + void SetGeoVertex(TVector3 vertex); + void SetVertexEnergy(Double_t vertexEner); + + AtFittedTrackOld &AddFittedTrack(std::unique_ptr ptr) + { + fFittedTrackArray.push_back(std::move(ptr)); + if (fFittedTrackArray.back()->GetTrackID() == -1) + fFittedTrackArray.back()->SetTrackID(fFittedTrackArray.size() - 1); + LOG(debug) << "Adding Track with ID " << fFittedTrackArray.back()->GetTrackID() << " to event " << fEventID; + + return *(fFittedTrackArray.back()); + } + + Double_t GetVertex(); + Double_t GetVertexEnergy(); + TVector3 GetGeoVertex(); + std::vector GetTrackArray(); + const FTrackVector &GetFittedTracks() const { return fFittedTrackArray; } + +private: + std::vector fTrackArray; + FTrackVector fFittedTrackArray; + Double_t fVertex{-10.0}; + Double_t fVertexEnergy{-10.0}; + TVector3 fGeoVertex; + + ClassDef(AtTrackingEventOld, 1); +}; + +#endif diff --git a/AtData/CMakeLists.txt b/AtData/CMakeLists.txt index 445cb359f..b3c0a5d1a 100644 --- a/AtData/CMakeLists.txt +++ b/AtData/CMakeLists.txt @@ -41,7 +41,6 @@ set(SRCS AtProtoQuadrant.cxx AtTrack.cxx AtContainerManip.cxx - AtMCResult.cxx AtFissionEvent.cxx @@ -55,7 +54,16 @@ set(SRCS AtPattern/AtPadPlaneElement.cxx AtFittedTrack.cxx - + AtFitResult.cxx + AtFitTrackResult.cxx + AtMCResult.cxx + + + # Deprecated things to be removed eventually. + AtTrackingEventOld.cxx + AtFittedTrackOld.cxx + AtMCResultOld.cxx + ) set(TEST_SRCS diff --git a/AtReconstruction/AtFitter/AtFitter.cxx b/AtReconstruction/AtFitter/AtFitter.cxx index bbb806185..9cca3032d 100644 --- a/AtReconstruction/AtFitter/AtFitter.cxx +++ b/AtReconstruction/AtFitter/AtFitter.cxx @@ -1,268 +1,11 @@ #include "AtFitter.h" -#include "AtHit.h" -#include "AtHitCluster.h" // for AtHitCluster -#include "AtTrack.h" - -#include // for Cartesian3D, operator-, PositionVector3D -#include // for DisplacementVector3D -#include - -#include -#include // for sqrt -#include // for operator<<, basic_ostream::operator<< -#include // for shared_ptr, __shared_ptr_access, __sha... -#include // for pair - -constexpr auto cRED = "\033[1;31m"; -constexpr auto cYELLOW = "\033[1;33m"; -constexpr auto cNORMAL = "\033[0m"; -constexpr auto cGREEN = "\033[1;32m"; - ClassImp(AtFITTER::AtFitter); -AtFITTER::AtFitter::AtFitter() = default; - -AtFITTER::AtFitter::~AtFitter() = default; - -std::tuple AtFITTER::AtFitter::GetMomFromBrho(Double_t M, Double_t Z, Double_t brho) -{ - - const Double_t M_Ener = M * 931.49401 / 1000.0; - Double_t p = brho * Z * (2.99792458 / 10.0); // In GeV - Double_t E = TMath::Sqrt(TMath::Power(p, 2) + TMath::Power(M_Ener, 2)) - M_Ener; - // std::cout << " Brho : " << brho << " - p : " << p << " - E : " << E << "\n"; - return std::make_tuple(p, E); -} - -Bool_t AtFITTER::AtFitter::FindVertexTrack(AtTrack *trA, AtTrack *trB) -{ - // Determination of first hit distance. NB: Assuming both tracks have the same angle sign - Double_t vertexA = 0.0; - Double_t vertexB = 0.0; - if (trA->GetGeoTheta() * TMath::RadToDeg() < 90) { - auto iniClusterA = trA->GetHitClusterArray()->back(); - auto iniClusterB = trB->GetHitClusterArray()->back(); - vertexA = 1000.0 - iniClusterA.GetPosition().Z(); - vertexB = 1000.0 - iniClusterB.GetPosition().Z(); - } else if (trA->GetGeoTheta() * TMath::RadToDeg() > 90) { - auto iniClusterA = trA->GetHitClusterArray()->front(); - auto iniClusterB = trB->GetHitClusterArray()->front(); - vertexA = iniClusterA.GetPosition().Z(); - vertexB = iniClusterB.GetPosition().Z(); - } - - return vertexA < vertexB; -} - -Bool_t AtFITTER::AtFitter::MergeTracks(std::vector *trackCandSource, std::vector *trackDest, - Bool_t enableSingleVertexTrack, Double_t clusterRadius, Double_t clusterDistance) +void AtFITTER::AtFitter::Reset() { - - Bool_t toMerge = kFALSE; - - Int_t addHitCnt = 0; - // Find the track closer to vertex - std::sort(trackCandSource->begin(), trackCandSource->end(), - [this](AtTrack *trA, AtTrack *trB) { return FindVertexTrack(trA, trB); }); - - // Track stitching from vertex - AtTrack *vertexTrack = *trackCandSource->begin(); - - if (enableSingleVertexTrack) { - - // Mark all tracks as merged - for (auto track : *trackCandSource) - track->SetIsMerged(kTRUE); - - trackDest->push_back(*vertexTrack); - return true; - } - - // Check if the candidate vertex track was merged - if (vertexTrack->GetIsMerged()) - return kFALSE; - else - vertexTrack->SetIsMerged(kTRUE); - - // If enabled, choose only the track closest to vertex (i.e. first one of the collection of candidates) - // TODO: Select by number of points - - for (auto it = trackCandSource->begin() + 1; it != trackCandSource->end(); ++it) { - // NB: These tracks were previously marked to merge. If merging fails they should be discarded. - AtTrack *trackToMerge = *(it); - toMerge = kFALSE; - - // Skip trackes flagged as merged - if (!trackToMerge->GetIsMerged()) { - trackToMerge->SetIsMerged(kTRUE); - } else - continue; - - Double_t endVertexZ = 0.0; - Double_t iniMergeZ = 0.0; - std::cout << " Trying to merge ... " - << "\n"; - std::cout << " Vertex track " << vertexTrack->GetTrackID() << " - Track to Merge " << trackToMerge->GetTrackID() - << "\n"; - // Check relative position between end and begin of each track using Hit Clusters - std::cout << " Vertex angle " << vertexTrack->GetGeoTheta() * TMath::RadToDeg() << "\n"; - if (vertexTrack->GetGeoTheta() * TMath::RadToDeg() < 90) { - auto endClusterVertex = vertexTrack->GetHitClusterArray()->front(); - auto iniClusterMerge = trackToMerge->GetHitClusterArray()->back(); - // Check separation and relative distance - endVertexZ = 1000.0 - endClusterVertex.GetPosition().Z(); - iniMergeZ = 1000.0 - iniClusterMerge.GetPosition().Z(); - - Double_t distance = std::sqrt((iniClusterMerge.GetPosition() - endClusterVertex.GetPosition()).Mag2()); - // std::cout << " Distance between tracks " << distance << "\n"; - // std::cout << " Ini Merge " << iniMergeZ << " - endVertexZ " << endVertexZ << "\n"; - if (((iniMergeZ + 10.0) > endVertexZ) && distance < 200) { - toMerge = kTRUE; - } - - } else if (vertexTrack->GetGeoTheta() * TMath::RadToDeg() > 90) { - auto endClusterVertex = vertexTrack->GetHitClusterArray()->back(); - auto iniClusterMerge = trackToMerge->GetHitClusterArray()->front(); - // Check separation and relative distance - endVertexZ = endClusterVertex.GetPosition().Z(); - iniMergeZ = iniClusterMerge.GetPosition().Z(); - Double_t distance = std::sqrt((iniClusterMerge.GetPosition() - endClusterVertex.GetPosition()).Mag2()); - // std::cout<<" Distance between tracks "< endVertexZ) && - distance < 100) { // NB: Distance between parts of the backward tracks is more critical - toMerge = kTRUE; - } - } - - if (toMerge) { - - std::cout << " --- Merging Succeeded! Vertex track " << vertexTrack->GetTrackID() << " - Track to Merge " - << trackToMerge->GetTrackID() << "\n"; - for (const auto &hit : trackToMerge->GetHitArray()) { - - vertexTrack->AddHit(hit->Clone()); // TODO: Look at code and see if this can be a move instead of a copy - ++addHitCnt; - } - - // Reclusterize after merging - vertexTrack->SortHitArrayTime(); - vertexTrack->ResetHitClusterArray(); - fTrackTransformer->ClusterizeSmooth3D( - *vertexTrack, clusterRadius, - clusterDistance); // NB: It can be removed if we force reclusterization for any track in the mina program - - // TODO: Check if phi recalculatio is needed - - } else { - std::cout << " --- Merging Failed ! Vertex track " << vertexTrack->GetTrackID() << " - Track to Merge " - << trackToMerge->GetTrackID() << "\n"; - } - } - - trackDest->push_back(*vertexTrack); - - return toMerge; -} -[[deprecated]] void -AtFITTER::AtFitter::MergeTracks(std::vector *trackCandSource, std::vector *trackJunkSource, - std::vector *trackDest, bool fitDirection, bool simulationConv) -{ - // DEPRECATED - // Track destination are the merged tracks. - // Track candidate source are the main tracks identified as candidates. - // Track junk source are the tracks from which clusters will be extracted. Boundary conditions are applied: Proximity - // in space, angle and center. - // NB: Only works for backward tracks - - Double_t trackDist = 20.0; // Distance between clusters in mm - Double_t angleSpread = 5.0; // Maximum angular spread between clusters - Double_t centerSpread = 20.0; // Maximim distance between track centers for merging. - Int_t minClusters = 3; - Int_t trackSize = 0; - - for (auto trackCand : *trackCandSource) { - Double_t thetaCand = trackCand.GetGeoTheta(); - auto &hitArrayCand = trackCand.GetHitArray(); - std::pair centerCand = trackCand.GetGeoCenter(); - - AtTrack track = trackCand; - Int_t jnkCnt = 0; - Int_t jnkHitCnt = 0; - - if (simulationConv) { - thetaCand = 180.0 - thetaCand * TMath::RadToDeg(); - - } else { - thetaCand = thetaCand * TMath::RadToDeg(); - } - - for (auto trackJunk : *trackJunkSource) { - Double_t thetaJunk = trackJunk.GetGeoTheta(); - auto &hitArrayJunk = trackJunk.GetHitArray(); - std::pair centerJunk = trackJunk.GetGeoCenter(); - - if (simulationConv) { - - thetaJunk = 180.0 - thetaJunk * TMath::RadToDeg(); - } else { - - thetaJunk = thetaJunk * TMath::RadToDeg(); - } - - if (thetaCand > 90) { - - /*if ((thetaCand + angleSpread) > thetaJunk && (thetaCand - angleSpread) < thetaJunk) { - for (auto hit : *hitArrayJunk) { - - track.AddHit(&hit); - ++jnkHitCnt; - } - }*/ - - } else if (thetaCand < 90 && thetaJunk < 90) { - - if ((thetaCand + angleSpread) > thetaJunk && (thetaCand - angleSpread) < thetaJunk) { // Check angle - std::cout << " Center cand : " << centerCand.first << " - " << centerCand.second << " " << thetaCand - << "\n"; - std::cout << " Center junk : " << centerJunk.first << " - " << centerJunk.second << " " << thetaJunk - << "\n"; - Double_t centerDistance = TMath::Sqrt(TMath::Power(centerCand.first - centerJunk.first, 2) + - TMath::Power(centerCand.second - centerJunk.second, 2)); - std::cout << " Distance " << centerDistance << "\n"; - if (centerDistance < 50) // Check quadrant - - for (const auto &hit : hitArrayJunk) { - - track.AddHit(hit->Clone()); - ++jnkHitCnt; - } - } - } - - ++jnkCnt; - - } // Junk track - - // track.SortClusterHitArrayZ(); - track.SortHitArrayTime(); - - // Prune if other tracks were added track - - if (trackJunkSource->size() > 0) { - Double_t pruneFraction = 10.0; - if (jnkHitCnt > hitArrayCand.size()) - pruneFraction = 4.0; // 25% - else - pruneFraction = 10.0; // 10% - - Int_t numHits = (Int_t)track.GetHitArray().size() / pruneFraction; - for (auto iHit = 0; iHit < numHits; ++iHit) - track.GetHitArray().pop_back(); - } - - trackDest->push_back(track); - - } // Source track + fPatternEvent = nullptr; + fFitResult = nullptr; + fRawEvent = nullptr; + fEvent = nullptr; } diff --git a/AtReconstruction/AtFitter/AtFitter.h b/AtReconstruction/AtFitter/AtFitter.h index 9d1e8274a..32233e28b 100644 --- a/AtReconstruction/AtFitter/AtFitter.h +++ b/AtReconstruction/AtFitter/AtFitter.h @@ -1,51 +1,71 @@ #ifndef AtFITTER_H #define AtFITTER_H -#include "AtTrackTransformer.h" +#include "AtFitResult.h" #include #include +#include #include +#include #include #include class AtDigiPar; class AtTrack; class AtFittedTrack; -class FairLogger; class TBuffer; class TClass; class TMemberInspector; - -namespace genfit { -class Track; -} // namespace genfit +class AtRawEvent; +class AtEvent; +class AtPatternEvent; namespace AtFITTER { class AtFitter : public TObject { +public: + using TrackResultPtr = std::unique_ptr; + using TrackResultsVector = std::vector; + using TrackResultsSet = + std::set>; + +protected: + // Pointer to the AtPatternEvent to be fitted. + AtPatternEvent *fPatternEvent{nullptr}; + + // Pointer to the AtFitResult object in which store the fit metadata. + AtFitResult *fFitResult{nullptr}; + + // Pointers to AtRawEvent and AtEvent. In case some specific fitter needs to access information + // in any of those branches. + AtRawEvent *fRawEvent{nullptr}; + AtEvent *fEvent{nullptr}; + + // Compare function that will be used to sort the fit results for a given track. + virtual bool CompareTrackFitsFunction(const TrackResultPtr &trackResultA, const TrackResultPtr &trackResultB) = 0; public: - AtFitter(); - virtual ~AtFitter(); - virtual std::vector> ProcessTracks(std::vector &tracks) = 0; + AtFitter() = default; + ~AtFitter() = default; + + virtual std::vector> ProcessEvent() = 0; virtual void Init() = 0; - void MergeTracks(std::vector *trackCandSource, std::vector *trackJunkSource, - std::vector *trackDest, bool fitDirection, bool simulationConv); - Bool_t MergeTracks(std::vector *trackCandSource, std::vector *trackDest, - Bool_t enableSingleVertexTrack, Double_t clusterRadius, Double_t clusterDistance); + // Mandatory to set. + void SetPatternEvent(AtPatternEvent *patternEvent) { fPatternEvent = patternEvent; } + void SetFitResult(AtFitResult *fitResult) { fFitResult = fitResult; } + + // Optional to set. + void SetRawEvent(AtRawEvent *rawEvent) { fRawEvent = rawEvent; } + void SetEvent(AtEvent *event) { fEvent = event; } + + // Reset pointers to AtPatternEvent, AtRawEvent and AtEvent. + void Reset(); protected: - FairLogger *fLogger{}; ///< logger pointer - AtDigiPar *fPar{}; ///< parameter container - std::unique_ptr fTrackTransformer{std::make_unique()}; - std::tuple - GetMomFromBrho(Double_t A, Double_t Z, - Double_t brho); ///< Returns momentum (in GeV) from Brho assuming M (amu) and Z; - Bool_t FindVertexTrack(AtTrack *trA, AtTrack *trB); ///< Lambda function to find track closer to vertex - ClassDef(AtFitter, 1); + ClassDef(AtFitter, 2); }; } // namespace AtFITTER diff --git a/AtReconstruction/AtFitter/AtFitterOld.cxx b/AtReconstruction/AtFitter/AtFitterOld.cxx new file mode 100644 index 000000000..a663bb92a --- /dev/null +++ b/AtReconstruction/AtFitter/AtFitterOld.cxx @@ -0,0 +1,269 @@ +#include "AtFitterOld.h" + +#include "AtHit.h" +#include "AtHitCluster.h" // for AtHitCluster +#include "AtTrack.h" + +#include // for Cartesian3D, operator-, PositionVector3D +#include // for DisplacementVector3D +#include + +#include +#include // for sqrt +#include // for operator<<, basic_ostream::operator<< +#include // for shared_ptr, __shared_ptr_access, __sha... +#include // for pair + +constexpr auto cRED = "\033[1;31m"; +constexpr auto cYELLOW = "\033[1;33m"; +constexpr auto cNORMAL = "\033[0m"; +constexpr auto cGREEN = "\033[1;32m"; + +ClassImp(AtFITTER::AtFitterOld); + +AtFITTER::AtFitterOld::AtFitterOld() = default; + +AtFITTER::AtFitterOld::~AtFitterOld() = default; + +std::tuple AtFITTER::AtFitterOld::GetMomFromBrho(Double_t M, Double_t Z, Double_t brho) +{ + + const Double_t M_Ener = M * 931.49401 / 1000.0; + Double_t p = brho * Z * (2.99792458 / 10.0); // In GeV + Double_t E = TMath::Sqrt(TMath::Power(p, 2) + TMath::Power(M_Ener, 2)) - M_Ener; + // std::cout << " Brho : " << brho << " - p : " << p << " - E : " << E << "\n"; + return std::make_tuple(p, E); +} + +Bool_t AtFITTER::AtFitterOld::FindVertexTrack(AtTrack *trA, AtTrack *trB) +{ + // Determination of first hit distance. NB: Assuming both tracks have the same angle sign + Double_t vertexA = 0.0; + Double_t vertexB = 0.0; + if (trA->GetGeoTheta() * TMath::RadToDeg() < 90) { + auto iniClusterA = trA->GetHitClusterArray()->back(); + auto iniClusterB = trB->GetHitClusterArray()->back(); + vertexA = 1000.0 - iniClusterA.GetPosition().Z(); + vertexB = 1000.0 - iniClusterB.GetPosition().Z(); + } else if (trA->GetGeoTheta() * TMath::RadToDeg() > 90) { + auto iniClusterA = trA->GetHitClusterArray()->front(); + auto iniClusterB = trB->GetHitClusterArray()->front(); + vertexA = iniClusterA.GetPosition().Z(); + vertexB = iniClusterB.GetPosition().Z(); + } + + return vertexA < vertexB; +} + +Bool_t AtFITTER::AtFitterOld::MergeTracks(std::vector *trackCandSource, std::vector *trackDest, + Bool_t enableSingleVertexTrack, Double_t clusterRadius, + Double_t clusterDistance) +{ + + Bool_t toMerge = kFALSE; + + Int_t addHitCnt = 0; + // Find the track closer to vertex + std::sort(trackCandSource->begin(), trackCandSource->end(), + [this](AtTrack *trA, AtTrack *trB) { return FindVertexTrack(trA, trB); }); + + // Track stitching from vertex + AtTrack *vertexTrack = *trackCandSource->begin(); + + if (enableSingleVertexTrack) { + + // Mark all tracks as merged + for (auto track : *trackCandSource) + track->SetIsMerged(kTRUE); + + trackDest->push_back(*vertexTrack); + return true; + } + + // Check if the candidate vertex track was merged + if (vertexTrack->GetIsMerged()) + return kFALSE; + else + vertexTrack->SetIsMerged(kTRUE); + + // If enabled, choose only the track closest to vertex (i.e. first one of the collection of candidates) + // TODO: Select by number of points + + for (auto it = trackCandSource->begin() + 1; it != trackCandSource->end(); ++it) { + // NB: These tracks were previously marked to merge. If merging fails they should be discarded. + AtTrack *trackToMerge = *(it); + toMerge = kFALSE; + + // Skip trackes flagged as merged + if (!trackToMerge->GetIsMerged()) { + trackToMerge->SetIsMerged(kTRUE); + } else + continue; + + Double_t endVertexZ = 0.0; + Double_t iniMergeZ = 0.0; + std::cout << " Trying to merge ... " + << "\n"; + std::cout << " Vertex track " << vertexTrack->GetTrackID() << " - Track to Merge " << trackToMerge->GetTrackID() + << "\n"; + // Check relative position between end and begin of each track using Hit Clusters + std::cout << " Vertex angle " << vertexTrack->GetGeoTheta() * TMath::RadToDeg() << "\n"; + if (vertexTrack->GetGeoTheta() * TMath::RadToDeg() < 90) { + auto endClusterVertex = vertexTrack->GetHitClusterArray()->front(); + auto iniClusterMerge = trackToMerge->GetHitClusterArray()->back(); + // Check separation and relative distance + endVertexZ = 1000.0 - endClusterVertex.GetPosition().Z(); + iniMergeZ = 1000.0 - iniClusterMerge.GetPosition().Z(); + + Double_t distance = std::sqrt((iniClusterMerge.GetPosition() - endClusterVertex.GetPosition()).Mag2()); + // std::cout << " Distance between tracks " << distance << "\n"; + // std::cout << " Ini Merge " << iniMergeZ << " - endVertexZ " << endVertexZ << "\n"; + if (((iniMergeZ + 10.0) > endVertexZ) && distance < 200) { + toMerge = kTRUE; + } + + } else if (vertexTrack->GetGeoTheta() * TMath::RadToDeg() > 90) { + auto endClusterVertex = vertexTrack->GetHitClusterArray()->back(); + auto iniClusterMerge = trackToMerge->GetHitClusterArray()->front(); + // Check separation and relative distance + endVertexZ = endClusterVertex.GetPosition().Z(); + iniMergeZ = iniClusterMerge.GetPosition().Z(); + Double_t distance = std::sqrt((iniClusterMerge.GetPosition() - endClusterVertex.GetPosition()).Mag2()); + // std::cout<<" Distance between tracks "< endVertexZ) && + distance < 100) { // NB: Distance between parts of the backward tracks is more critical + toMerge = kTRUE; + } + } + + if (toMerge) { + + std::cout << " --- Merging Succeeded! Vertex track " << vertexTrack->GetTrackID() << " - Track to Merge " + << trackToMerge->GetTrackID() << "\n"; + for (const auto &hit : trackToMerge->GetHitArray()) { + + vertexTrack->AddHit(hit->Clone()); // TODO: Look at code and see if this can be a move instead of a copy + ++addHitCnt; + } + + // Reclusterize after merging + vertexTrack->SortHitArrayTime(); + vertexTrack->ResetHitClusterArray(); + fTrackTransformer->ClusterizeSmooth3D( + *vertexTrack, clusterRadius, + clusterDistance); // NB: It can be removed if we force reclusterization for any track in the mina program + + // TODO: Check if phi recalculatio is needed + + } else { + std::cout << " --- Merging Failed ! Vertex track " << vertexTrack->GetTrackID() << " - Track to Merge " + << trackToMerge->GetTrackID() << "\n"; + } + } + + trackDest->push_back(*vertexTrack); + + return toMerge; +} +[[deprecated]] void +AtFITTER::AtFitterOld::MergeTracks(std::vector *trackCandSource, std::vector *trackJunkSource, + std::vector *trackDest, bool fitDirection, bool simulationConv) +{ + // DEPRECATED + // Track destination are the merged tracks. + // Track candidate source are the main tracks identified as candidates. + // Track junk source are the tracks from which clusters will be extracted. Boundary conditions are applied: Proximity + // in space, angle and center. + // NB: Only works for backward tracks + + Double_t trackDist = 20.0; // Distance between clusters in mm + Double_t angleSpread = 5.0; // Maximum angular spread between clusters + Double_t centerSpread = 20.0; // Maximim distance between track centers for merging. + Int_t minClusters = 3; + Int_t trackSize = 0; + + for (auto trackCand : *trackCandSource) { + Double_t thetaCand = trackCand.GetGeoTheta(); + auto &hitArrayCand = trackCand.GetHitArray(); + std::pair centerCand = trackCand.GetGeoCenter(); + + AtTrack track = trackCand; + Int_t jnkCnt = 0; + Int_t jnkHitCnt = 0; + + if (simulationConv) { + thetaCand = 180.0 - thetaCand * TMath::RadToDeg(); + + } else { + thetaCand = thetaCand * TMath::RadToDeg(); + } + + for (auto trackJunk : *trackJunkSource) { + Double_t thetaJunk = trackJunk.GetGeoTheta(); + auto &hitArrayJunk = trackJunk.GetHitArray(); + std::pair centerJunk = trackJunk.GetGeoCenter(); + + if (simulationConv) { + + thetaJunk = 180.0 - thetaJunk * TMath::RadToDeg(); + } else { + + thetaJunk = thetaJunk * TMath::RadToDeg(); + } + + if (thetaCand > 90) { + + /*if ((thetaCand + angleSpread) > thetaJunk && (thetaCand - angleSpread) < thetaJunk) { + for (auto hit : *hitArrayJunk) { + + track.AddHit(&hit); + ++jnkHitCnt; + } + }*/ + + } else if (thetaCand < 90 && thetaJunk < 90) { + + if ((thetaCand + angleSpread) > thetaJunk && (thetaCand - angleSpread) < thetaJunk) { // Check angle + std::cout << " Center cand : " << centerCand.first << " - " << centerCand.second << " " << thetaCand + << "\n"; + std::cout << " Center junk : " << centerJunk.first << " - " << centerJunk.second << " " << thetaJunk + << "\n"; + Double_t centerDistance = TMath::Sqrt(TMath::Power(centerCand.first - centerJunk.first, 2) + + TMath::Power(centerCand.second - centerJunk.second, 2)); + std::cout << " Distance " << centerDistance << "\n"; + if (centerDistance < 50) // Check quadrant + + for (const auto &hit : hitArrayJunk) { + + track.AddHit(hit->Clone()); + ++jnkHitCnt; + } + } + } + + ++jnkCnt; + + } // Junk track + + // track.SortClusterHitArrayZ(); + track.SortHitArrayTime(); + + // Prune if other tracks were added track + + if (trackJunkSource->size() > 0) { + Double_t pruneFraction = 10.0; + if (jnkHitCnt > hitArrayCand.size()) + pruneFraction = 4.0; // 25% + else + pruneFraction = 10.0; // 10% + + Int_t numHits = (Int_t)track.GetHitArray().size() / pruneFraction; + for (auto iHit = 0; iHit < numHits; ++iHit) + track.GetHitArray().pop_back(); + } + + trackDest->push_back(track); + + } // Source track +} diff --git a/AtReconstruction/AtFitter/AtFitterOld.h b/AtReconstruction/AtFitter/AtFitterOld.h new file mode 100644 index 000000000..1f08e1252 --- /dev/null +++ b/AtReconstruction/AtFitter/AtFitterOld.h @@ -0,0 +1,53 @@ +#ifndef AtFITTEROLD_H +#define AtFITTEROLD_H + +#include "AtTrackTransformer.h" + +#include +#include + +#include +#include +#include + +class AtDigiPar; +class AtTrack; +class AtFittedTrack; +class FairLogger; +class TBuffer; +class TClass; +class TMemberInspector; + +namespace genfit { +class Track; +} // namespace genfit + +namespace AtFITTER { + +class [[deprecated]] AtFitterOld : public TObject { + +public: + AtFitterOld(); + virtual ~AtFitterOld(); + virtual std::vector> ProcessTracks(std::vector &tracks) = 0; + virtual void Init() = 0; + + void MergeTracks(std::vector *trackCandSource, std::vector *trackJunkSource, + std::vector *trackDest, bool fitDirection, bool simulationConv); + Bool_t MergeTracks(std::vector *trackCandSource, std::vector *trackDest, + Bool_t enableSingleVertexTrack, Double_t clusterRadius, Double_t clusterDistance); + +protected: + FairLogger *fLogger{}; ///< logger pointer + AtDigiPar *fPar{}; ///< parameter container + std::unique_ptr fTrackTransformer{std::make_unique()}; + std::tuple + GetMomFromBrho(Double_t A, Double_t Z, + Double_t brho); ///< Returns momentum (in GeV) from Brho assuming M (amu) and Z; + Bool_t FindVertexTrack(AtTrack *trA, AtTrack *trB); ///< Lambda function to find track closer to vertex + ClassDef(AtFitterOld, 1); +}; + +} // namespace AtFITTER + +#endif diff --git a/AtReconstruction/AtFitter/AtMCFitter.cxx b/AtReconstruction/AtFitter/AtMCFitter.cxx index bef820f99..ea455eccd 100644 --- a/AtReconstruction/AtFitter/AtMCFitter.cxx +++ b/AtReconstruction/AtFitter/AtMCFitter.cxx @@ -30,7 +30,7 @@ namespace MCFitter { AtMCFitter::AtMCFitter(SimPtr sim, ClusterPtr cluster, PulsePtr pulse) : fMap(pulse->GetMap()), fSim(move(sim)), fClusterize(move(cluster)), fPulse(move(pulse)), - fResults([](const AtMCResult &a, const AtMCResult &b) { return a.fObjective < b.fObjective; }) + fResults([](const AtMCResult &a, const AtMCResult &b) { return a.GetChi2() < b.GetChi2(); }) { } @@ -84,7 +84,7 @@ void AtMCFitter::RunIterRange(int startIter, int numIter, AtPulse *pulse) double obj = ObjectiveFunction(*fCurrentEvent, idx, result); result.fIterNum = idx; - result.fObjective = obj; + result.SetChi2(obj); // result.Print(); { std::lock_guard lk(fResultMutex); diff --git a/AtReconstruction/AtFitter/AtMCFitterOld.cxx b/AtReconstruction/AtFitter/AtMCFitterOld.cxx new file mode 100644 index 000000000..87aa7985f --- /dev/null +++ b/AtReconstruction/AtFitter/AtMCFitterOld.cxx @@ -0,0 +1,216 @@ +#include "AtMCFitterOld.h" +// IWYU pragma: no_include + +#include "AtClusterize.h" // for AtClusterize +#include "AtDigiPar.h" // for AtDigiPar +#include "AtEvent.h" // for AtEvent +#include "AtMCResultOld.h" +#include "AtPSA.h" // for AtPSA +#include "AtParameterDistribution.h" +#include "AtPatternEvent.h" // for AtPatternEvent +#include "AtPulse.h" // for AtPulse +#include "AtRawEvent.h" // for AtRawEvent +#include "AtSimpleSimulation.h" // for AtSimpleSimulation +#include "AtSimulatedPoint.h" // IWYU pragma: keep +#include "AtSpaceChargeModel.h" + +#include // for LOG, Logger +#include // for FairParSet +#include // for FairRunAna +#include // for FairRuntimeDb + +#include + +#include // for max +#include +#include +#include +using std::move; +namespace MCFitter { + +AtMCFitterOld::AtMCFitterOld(SimPtr sim, ClusterPtr cluster, PulsePtr pulse) + : fMap(pulse->GetMap()), fSim(move(sim)), fClusterize(move(cluster)), fPulse(move(pulse)), + fResults([](const AtMCResultOld &a, const AtMCResultOld &b) { return a.fObjective < b.fObjective; }) +{ +} + +AtMCFitterOld::ParamPtr AtMCFitterOld::GetParameter(const std::string &name) const +{ + if (fParameters.find(name) != fParameters.end()) { + return fParameters.at(name); + } + return nullptr; +} + +void AtMCFitterOld::SetNumThreads(int num) +{ + if (num > 1) + ROOT::EnableThreadSafety(); + fNumThreads = num; +} + +void AtMCFitterOld::Init() +{ + CreateParamDistros(); + + FairRunAna *ana = FairRunAna::Instance(); + FairRuntimeDb *rtdb = ana->GetRuntimeDb(); + fPar = dynamic_cast(rtdb->getContainer("AtDigiPar")); + + fPulse->SetParameters(fPar); + fClusterize->GetParameters(fPar); + if (fPSA) + fPSA->Init(); + if (fSim->GetSpaceChargeModel()) + fSim->GetSpaceChargeModel()->LoadParameters(fPar); + + fThPulse.resize(fNumThreads); + for (int i = 0; i < fNumThreads; ++i) + fThPulse[i] = fPulse->Clone(); +} + +void AtMCFitterOld::RunIterRange(int startIter, int numIter, AtPulse *pulse) +{ + // Here we should copy each thread their own version of the clusterize, pulse, and simulation + // objects (only if the number of threads is greater than 1). Needs to be deep copies + + for (int i = 0; i < numIter; ++i) { + + int idx = startIter + i; + auto result = DefineEvent(); + auto mcPoints = SimulateEvent(result); + + DigitizeEvent(mcPoints, idx, pulse); + double obj = ObjectiveFunction(*fCurrentEvent, idx, result); + + result.fIterNum = idx; + result.fObjective = obj; + // result.Print(); + { + std::lock_guard lk(fResultMutex); + fResults.insert(result); + } + } + LOG(debug) << "Done with run iter range"; +} + +void AtMCFitterOld::Exec(const AtPatternEvent &event) +{ + fRawEventArray.clear(); + fEventArray.clear(); + fResults.clear(); + + SetParamDistributions(event); + + // Set the conditions for simulating the event + fCurrentEvent = &event; + + // Make sure the event arrays are large enough so no resizing will happen + fRawEventArray.resize(fNumIter); + fEventArray.resize(fNumIter); + + for (int i = 0; i < fNumRounds; ++i) { + RunRound(); + RecenterParamDistributions(); + } +} +void AtMCFitterOld::RunRound() +{ + // Begining of round + auto start = std::chrono::high_resolution_clock::now(); + + // Get what iterations to do on what thread. + std::vector> threadParam; + int iterPerTh = fNumIter / fNumThreads; + for (int i = 0; i < fNumThreads; ++i) + threadParam.emplace_back(0, iterPerTh); + for (int i = 0; i < fNumIter % fNumThreads; ++i) + threadParam[i].second++; + for (int i = 1; i < fNumThreads; ++i) + threadParam[i].first = threadParam[i - 1].first + threadParam[i - 1].second; + + for (int i = 0; i < threadParam.size(); ++i) { + LOG(info) << i << ": " << threadParam[i].first << " " << threadParam[i].second; + } + + std::vector threads; + for (int i = 0; i < fNumThreads; ++i) { + LOG(debug) << "Creating thread " << i << " with " << threadParam[i].first << " " << threadParam[i].second + << " and " << fPulse.get(); + + // Spawn a thread to call RunIterRange. + threads.emplace_back( + [this](std::pair param, AtPulse *pulse) { this->RunIterRange(param.first, param.second, pulse); }, + threadParam[i], fThPulse[i].get()); + } + + // Wait for all threads to finish + for (auto &th : threads) + th.join(); + + auto stop = std::chrono::high_resolution_clock::now(); + + if (fTimeEvent) + LOG(info) << "Simulation of " << fNumIter << " events took " + << std::chrono::duration_cast(stop - start).count() << " ms."; +} + +int AtMCFitterOld::DigitizeEvent(const TClonesArray &points, int idx, AtPulse *pulse) +{ + // Event has been simulated and is sitting in the fSim + auto vec = fClusterize->ProcessEvent(points); + LOG(debug) << "Digitizing event at " << idx; + + fRawEventArray[idx] = pulse->GenerateEvent(vec); + + if (fPSA) { + LOG(debug) << "Running PSA at " << idx; + fEventArray[idx] = fPSA->Analyze(fRawEventArray[idx]); + } + LOG(debug) << "Done digitizing event at " << idx; + return idx; +} + +/** + * Fill the TClonesArray in order of smallest to largest chi2. + */ +void AtMCFitterOld::FillResultArrays(TClonesArray &resultArray, TClonesArray &simEvent, TClonesArray &simRawEvent) +{ + resultArray.Delete(); + simEvent.Delete(); + simRawEvent.Delete(); + + for (auto &res : fResults) { + + int clonesIdx = resultArray.GetEntries(); + int eventIdx = res.fIterNum; + LOG(debug) << "Filling iteration " << eventIdx << " at index " << resultArray.GetEntries(); + + new (resultArray[clonesIdx]) AtMCResultOld(std::move(res)); + if (clonesIdx < fNumEventsToSave) { + new (simEvent[clonesIdx]) AtEvent(std::move(fEventArray[eventIdx])); + new (simRawEvent[clonesIdx]) AtRawEvent(std::move(fRawEventArray[eventIdx])); + } + } + + fEventArray.clear(); + fRawEventArray.clear(); +} + +AtMCResultOld AtMCFitterOld::DefineEvent() +{ + AtMCResultOld result; + for (auto &[name, distro] : fParameters) + result.fParameters[name] = distro->Sample(); + return result; +} +void AtMCFitterOld::RecenterParamDistributions() +{ + for (auto &[name, distro] : fParameters) { + AtMCResultOld result = *fResults.begin(); + distro->SetMean(result.fParameters[name]); + distro->TruncateSpace(); + } +} + +} // namespace MCFitter diff --git a/AtReconstruction/AtFitter/AtMCFitterOld.h b/AtReconstruction/AtFitter/AtMCFitterOld.h new file mode 100644 index 000000000..9cbdd6995 --- /dev/null +++ b/AtReconstruction/AtFitter/AtMCFitterOld.h @@ -0,0 +1,136 @@ +#ifndef ATMCFITTEROLD_H +#define ATMCFITTEROLD_H + +#include "AtEvent.h" +#include "AtMCResultOld.h" // for AtMCResult +#include "AtRawEvent.h" + +#include // for TClonesArray + +#include // for function +#include // for map +#include // for shared_ptr +#include // for mutex +#include // for set +#include // for string +#include // for pair +#include // for vector + +class AtBaseEvent; // lines 13-13 +class AtClusterize; // lines 15-15 +class AtMap; // lines 18-18 +class AtPatternEvent; // lines 12-12 +class AtPulse; // lines 16-16 +class AtSimpleSimulation; // lines 14-14 +class AtDigiPar; +class AtPSA; + +namespace MCFitter { +class AtParameterDistribution; + +class [[deprecated]] AtMCFitterOld { +protected: + using ParamPtr = std::shared_ptr; + using SimPtr = std::shared_ptr; + + using ClusterPtr = std::shared_ptr; + using PulsePtr = std::shared_ptr; + using MapPtr = std::shared_ptr; + using PsaPtr = std::shared_ptr; + using ObjPair = std::pair; //< Iteration number and objective function value + + std::map fParameters; + + MapPtr fMap; + SimPtr fSim; + ClusterPtr fClusterize; + PulsePtr fPulse; + PsaPtr fPSA{nullptr}; + + int fNumIter{1}; + int fNumRounds{1}; + int fNumEventsToSave{10}; + bool fTimeEvent{false}; + int fNumThreads{1}; + + // Things used by threads excecuting that are either expensive to create and delete + // or unaccessable due to FairRoot design choices + const AtPatternEvent *fCurrentEvent{nullptr}; + std::vector fThPulse; //< Cached because it is expensive to create and delete. + const AtDigiPar *fPar{nullptr}; // fRawEventArray; + std::vector fEventArray; + + /** Things below here need to be written to by threads and will be locked using a shared mutex ***/ + /// Store the iteration number sorted by lowest objective funtion + std::mutex fResultMutex; + std::set> fResults; + +public: + AtMCFitterOld(SimPtr sim, ClusterPtr cluster, PulsePtr pulse); + virtual ~AtMCFitterOld() = default; + + void Init(); + void SetPSA(PsaPtr psa) { fPSA = psa; } + void Exec(const AtPatternEvent &event); + + ParamPtr GetParameter(const std::string &name) const; + void FillResultArrays(TClonesArray &resultArray, TClonesArray &simEvent, TClonesArray &simRawEvent); + void SetNumIter(int iter) { fNumIter = iter; } + + /// Set number of times to run fNumIter iterations and then re-center and truncate the parameter space. + void SetNumRounds(int rounds) { fNumRounds = rounds; } + void SetTimeEvent(bool val) { fTimeEvent = val; } + void SetNumEventsToSave(int num) { fNumEventsToSave = num; } + void SetNumThreads(int num); + +protected: + void RunRound(); + void RunIterRange(int startIter, int numIter, AtPulse *pulse); + + /** + *@brief Create the parameter distributions to use for the fit. + */ + virtual void CreateParamDistros() = 0; + + /** + * @brief Set parameter distributions (mean/spread) from the event. + */ + virtual void SetParamDistributions(const AtPatternEvent &event) = 0; + + /** + * @brief This is the thing we are minimizing between events (SimEventID is index in TClonesArray) + */ + virtual double ObjectiveFunction(const AtBaseEvent &expEvent, int SimEventID, AtMCResultOld &definition) = 0; + + /** + * Simulate an event using the parameters in the passed AtMCResult class and return an array of + * the AtMCPoints to then digitize. + */ + virtual TClonesArray SimulateEvent(AtMCResultOld &definition) = 0; + + /** + * Sample parameter distributions and constrain the system to simulate an event. + * The parameters in AtMCResult will be used to then simulate an event. + * This function calls Sample() on all the parameter distributions and saves them. + */ + virtual AtMCResultOld DefineEvent(); + + /** + * Recenter the parameter distributions around the best result and truncate the parameter space. + */ + virtual void RecenterParamDistributions(); + + /** + * Create the AtRawEvent and AtEvent from fSim + * returns the index of the event in the TClonesArray + */ + int DigitizeEvent(const TClonesArray &points, int idx, AtPulse *pulse); +}; + +} // namespace MCFitter + +#endif // ATMCFITTEROLD_H diff --git a/AtReconstruction/AtFitter/AtMCFitterTaskOld.cxx b/AtReconstruction/AtFitter/AtMCFitterTaskOld.cxx new file mode 100644 index 000000000..10df98724 --- /dev/null +++ b/AtReconstruction/AtFitter/AtMCFitterTaskOld.cxx @@ -0,0 +1,51 @@ +#include "AtMCFitterTaskOld.h" + +#include "AtMCFitterOld.h" +#include "AtPatternEvent.h" + +#include // for LOG, Logger +#include // for FairRootManager + +#include // for TClonesArray +#include + +#include + +AtMCFitterTaskOld::AtMCFitterTaskOld(std::shared_ptr fitter) + : fFitter(std::move(fitter)), fResultArray("MCFitter::AtMCResultOld"), fSimEventArray("AtEvent"), + fSimRawEventArray("AtRawEvent") +{ +} + +InitStatus AtMCFitterTaskOld::Init() +{ + LOG(debug) << "Initialing fitter"; + fFitter->Init(); + + FairRootManager *ioman = FairRootManager::Instance(); + ioman->Register("SimEvent", "cbmsim", &fSimEventArray, fSaveEvent); + ioman->Register("SimRawEvent", "cbmsim", &fSimRawEventArray, fSaveRawEvent); + ioman->Register("AtMCResultOld", "cbmsim", &fResultArray, fSaveResult); + + fPatternArray = dynamic_cast(ioman->GetObject(fPatternBranchName)); + if (fPatternArray == nullptr) + LOG(fatal) << "Failed to load branch " << fPatternBranchName; + + LOG(debug) << "Done with sim init"; + return kSUCCESS; +} + +void AtMCFitterTaskOld::Exec(Option_t *) +{ + LOG(debug) << "Exec"; + auto patEvent = dynamic_cast(fPatternArray->At(0)); + if (!patEvent->IsGood()) + return; + + fFitter->Exec(*patEvent); + fResultArray.Delete(); + fSimEventArray.Delete(); + fSimRawEventArray.Delete(); + + fFitter->FillResultArrays(fResultArray, fSimEventArray, fSimRawEventArray); +} diff --git a/AtReconstruction/AtFitter/AtMCFitterTaskOld.h b/AtReconstruction/AtFitter/AtMCFitterTaskOld.h new file mode 100644 index 000000000..29d115203 --- /dev/null +++ b/AtReconstruction/AtFitter/AtMCFitterTaskOld.h @@ -0,0 +1,43 @@ +#ifndef ATMCFITTERTASKOLD_H +#define ATMCFITTERTASKOLD_H + +#include + +#include // for Option_t +#include +#include // for TString + +#include // for shared_ptr + +namespace MCFitter { +class AtMCFitterOld; +} + +class [[deprecated]] AtMCFitterTaskOld : public FairTask { + + std::shared_ptr fFitter; //! + TString fPatternBranchName{"AtPatternEvent"}; + TClonesArray *fPatternArray{nullptr}; + + TClonesArray fResultArray; //< Output of task + TClonesArray fSimEventArray; + TClonesArray fSimRawEventArray; + + Bool_t fSaveResult{true}; + Bool_t fSaveEvent{false}; + Bool_t fSaveRawEvent{false}; + +public: + AtMCFitterTaskOld(std::shared_ptr fitter); + + InitStatus Init() override; + void Exec(Option_t *option = "") override; + void Finish() override{}; + + void SetPatternBranchName(TString name) { fPatternBranchName = name; } + void SetSaveResult(bool val) { fSaveResult = val; } + void SetSaveEvent(bool val) { fSaveEvent = val; } + void SetSaveRawEvent(bool val) { fSaveRawEvent = val; } +}; + +#endif // ATMCFITTERTASKOLD_H diff --git a/AtReconstruction/AtFitterTask.cxx b/AtReconstruction/AtFitterTask.cxx index 8bef7da49..5f73e967f 100644 --- a/AtReconstruction/AtFitterTask.cxx +++ b/AtReconstruction/AtFitterTask.cxx @@ -1,10 +1,11 @@ #include "AtFitterTask.h" #include "AtDigiPar.h" +#include "AtEvent.h" #include "AtFitter.h" -#include "AtGenfit.h" #include "AtParsers.h" #include "AtPatternEvent.h" +#include "AtRawEvent.h" #include "AtTrackingEvent.h" #include @@ -15,7 +16,6 @@ #include #include -#include #include #include @@ -27,7 +27,8 @@ ClassImp(AtFitterTask); AtFitterTask::AtFitterTask(std::unique_ptr fitter) : fInputBranchName("AtPatternEvent"), fOutputBranchName("AtTrackingEvent"), fIsPersistence(kFALSE), - fTrackingEventArray(TClonesArray("AtTrackingEvent", 1)), fFitter(std::move(fitter)) + fTrackingEventArray(TClonesArray("AtTrackingEvent", 1)), fFitter(std::move(fitter)), fRawEventBranchName(""), + fEventBranchName(""), fFitResultBranchName("") { } @@ -46,6 +47,22 @@ void AtFitterTask::SetOutputBranch(TString branchName) fOutputBranchName = branchName; } +void AtFitterTask::SetRawEventBranch(TString branchName) +{ + fRawEventBranchName = branchName; +} + +void AtFitterTask::SetEventBranch(TString branchName) +{ + fEventBranchName = branchName; +} + +void AtFitterTask::SetFitResultBranch(TString branchName) +{ + fFitResultBranchName = branchName; + fSaveFitResult = true; +} + InitStatus AtFitterTask::Init() { FairRootManager *ioMan = FairRootManager::Instance(); @@ -54,13 +71,24 @@ InitStatus AtFitterTask::Init() return kERROR; } - fPatternEventArray = dynamic_cast(ioMan->GetObject("AtPatternEvent")); + fPatternEventArray = dynamic_cast(ioMan->GetObject(fInputBranchName)); if (fPatternEventArray == nullptr) { LOG(error) << "Cannot find AtPatternEvent array!"; return kERROR; } ioMan->Register(fOutputBranchName, "AtTPC", &fTrackingEventArray, fIsPersistence); + ioMan->Register(fFitResultBranchName, "AtTPC", &fFitResultArray, fIsPersistence && fSaveFitResult); + + fRawEventArray = dynamic_cast(ioMan->GetObject(fRawEventBranchName)); + if (fRawEventArray == nullptr) { + LOG(info) << "AtRawEvent branch name was not set. No AtRawEvent will be passed to the fitter."; + } + + fEventArray = dynamic_cast(ioMan->GetObject(fEventBranchName)); + if (fEventArray == nullptr) { + LOG(info) << "AtEvent branch name was not set. No AtEvent will be passed to the fitter."; + } return kSUCCESS; } @@ -84,22 +112,39 @@ void AtFitterTask::SetParContainers() void AtFitterTask::Exec(Option_t *option) { + fFitter->Reset(); + if (fPatternEventArray->GetEntriesFast() == 0) return; + // If there is AtRawEvent available, pass it to the fitter. + if (fRawEventArray) { + AtRawEvent *rawEvent = dynamic_cast(fRawEventArray->At(0)); + fFitter->SetRawEvent(rawEvent); + } + + // If there is AtEvent available, pass it to the fitter. + if (fEventArray) { + AtEvent *event = dynamic_cast(fEventArray->At(0)); + fFitter->SetEvent(event); + } + fTrackingEventArray.Delete(); auto trackingEvent = dynamic_cast(fTrackingEventArray.ConstructedAt(0)); + auto fitResult = dynamic_cast(fFitResultArray.ConstructedAt(0)); - std::cout << " Event Counter " << fEventCnt << "\n"; + LOG(info) << " AtFitterTask::Exec() : Fitting event " << fEventCnt; - AtPatternEvent &patternEvent = *(dynamic_cast(fPatternEventArray->At(0))); - std::vector &tracks = patternEvent.GetTrackCand(); - std::cout << " AtFitterTask:Exec - Number of candidate tracks : " << tracks.size() << "\n"; + AtPatternEvent *patternEvent = dynamic_cast(fPatternEventArray->At(0)); + std::vector &tracks = patternEvent->GetTrackCand(); + LOG(info) << " AtFitterTask::Exec() : Number of candidate tracks : " << tracks.size(); - auto fittedTracks = fFitter->ProcessTracks(tracks); + fFitter->SetPatternEvent(patternEvent); + fFitter->SetFitResult(fitResult); + auto fittedTracks = fFitter->ProcessEvent(); - std::cout << " Number of fitted tracks " << fittedTracks.size() << "\n"; + LOG(info) << " AtFitterTask::Exec() : Number of fitted tracks : " << fittedTracks.size(); for (auto &fittedTrack : fittedTracks) trackingEvent->AddFittedTrack(std::move(fittedTrack)); diff --git a/AtReconstruction/AtFitterTask.h b/AtReconstruction/AtFitterTask.h index dc3256b12..74a706227 100644 --- a/AtReconstruction/AtFitterTask.h +++ b/AtReconstruction/AtFitterTask.h @@ -15,9 +15,8 @@ #include #include +#include -#include "EventDisplay.h" -#include "Exception.h" #include "FairLogger.h" #include "FairRootManager.h" #include "FairRun.h" @@ -41,21 +40,22 @@ class AtTrackTransformer; namespace AtFITTER { class AtFitter; } // namespace AtFITTER -namespace genfit { -class Track; -} // namespace genfit class AtFitterTask : public FairTask { public: - // AtFitterTask(); - ~AtFitterTask() = default; AtFitterTask(std::unique_ptr fitter); + ~AtFitterTask() = default; void SetInputBranch(TString branchName); void SetOutputBranch(TString branchName); void SetPersistence(Bool_t value = kTRUE); + void SetRawEventBranch(TString branchName); + void SetEventBranch(TString branchName); + + void SetFitResultBranch(TString branchName); + virtual InitStatus Init(); virtual void SetParContainers(); virtual void Exec(Option_t *opt); @@ -73,7 +73,18 @@ class AtFitterTask : public FairTask { std::size_t fEventCnt{0}; - ClassDef(AtFitterTask, 1); + // Include the option to input AtRawEvent and AtEvent in case some specific AtFitter needs it. + TString fRawEventBranchName; + TString fEventBranchName; + TClonesArray *fRawEventArray; + TClonesArray *fEventArray; + + // Include the option to store all the fit metadata in a AtFitResult branch. + TString fFitResultBranchName; + TClonesArray fFitResultArray; + bool fSaveFitResult{false}; + + ClassDef(AtFitterTask, 2); }; #endif diff --git a/AtReconstruction/AtFitterTaskOld.cxx b/AtReconstruction/AtFitterTaskOld.cxx new file mode 100644 index 000000000..7c5be94e0 --- /dev/null +++ b/AtReconstruction/AtFitterTaskOld.cxx @@ -0,0 +1,108 @@ +#include "AtFitterTaskOld.h" + +#include "AtDigiPar.h" +#include "AtFitter.h" +#include "AtGenfit.h" +#include "AtParsers.h" +#include "AtPatternEvent.h" +#include "AtTrackingEvent.h" + +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +class AtTrack; +class AtFittedTrackOld; + +ClassImp(AtFitterTaskOld); + +AtFitterTaskOld::AtFitterTaskOld(std::unique_ptr fitter) + : fInputBranchName("AtPatternEvent"), fOutputBranchName("AtTrackingEvent"), fIsPersistence(kFALSE), + fTrackingEventArray(TClonesArray("AtTrackingEvent", 1)), fFitter(std::move(fitter)) +{ +} + +void AtFitterTaskOld::SetPersistence(Bool_t value) +{ + fIsPersistence = value; +} + +void AtFitterTaskOld::SetInputBranch(TString branchName) +{ + fInputBranchName = branchName; +} + +void AtFitterTaskOld::SetOutputBranch(TString branchName) +{ + fOutputBranchName = branchName; +} + +InitStatus AtFitterTaskOld::Init() +{ + FairRootManager *ioMan = FairRootManager::Instance(); + if (ioMan == nullptr) { + LOG(error) << "Cannot find RootManager!"; + return kERROR; + } + + fPatternEventArray = dynamic_cast(ioMan->GetObject("AtPatternEvent")); + if (fPatternEventArray == nullptr) { + LOG(error) << "Cannot find AtPatternEvent array!"; + return kERROR; + } + + ioMan->Register(fOutputBranchName, "AtTPC", &fTrackingEventArray, fIsPersistence); + + return kSUCCESS; +} + +void AtFitterTaskOld::SetParContainers() +{ + LOG(debug) << "SetParContainers of AtFitterTask"; + + FairRun *run = FairRun::Instance(); + if (!run) + LOG(fatal) << "No analysis run!"; + + FairRuntimeDb *db = run->GetRuntimeDb(); // NOLINT + if (!db) + LOG(fatal) << "No runtime database!"; + + fPar = (AtDigiPar *)db->getContainer("AtDigiPar"); // NOLINT + if (!fPar) + LOG(fatal) << "AtDigiPar not found!!"; +} + +void AtFitterTaskOld::Exec(Option_t *option) +{ + if (fPatternEventArray->GetEntriesFast() == 0) + return; + + fTrackingEventArray.Delete(); + + auto trackingEvent = dynamic_cast(fTrackingEventArray.ConstructedAt(0)); + + std::cout << " Event Counter " << fEventCnt << "\n"; + + AtPatternEvent &patternEvent = *(dynamic_cast(fPatternEventArray->At(0))); + std::vector &tracks = patternEvent.GetTrackCand(); + std::cout << " AtFitterTask:Exec - Number of candidate tracks : " << tracks.size() << "\n"; + + auto fittedTracks = fFitter->ProcessTracks(tracks); + + std::cout << " Number of fitted tracks " << fittedTracks.size() << "\n"; + + for (auto &fittedTrack : fittedTracks) + trackingEvent->AddFittedTrack(std::move(fittedTrack)); + + ++fEventCnt; +} diff --git a/AtReconstruction/AtFitterTaskOld.h b/AtReconstruction/AtFitterTaskOld.h new file mode 100644 index 000000000..c627a37b9 --- /dev/null +++ b/AtReconstruction/AtFitterTaskOld.h @@ -0,0 +1,79 @@ +/********************************************************************* + * Fitter Task AtFitterTask.hh * + * Author: Y. Ayyad ayyadlim@frib.msu.edu * + * Log: 3/10/2021 * + * * + *********************************************************************/ + +#ifndef ATFITTERTASKOLD +#define ATFITTERTASKOLD + +#include "AtFormat.h" +#include "AtKinematics.h" +#include "AtParsers.h" + +#include + +#include + +#include "EventDisplay.h" +#include "Exception.h" +#include "FairLogger.h" +#include "FairRootManager.h" +#include "FairRun.h" +#include "FairRunAna.h" + +#include +#include +#include + +class AtDigiPar; +class FairLogger; +class TBuffer; +class TClass; +class TClonesArray; +class TMemberInspector; +class AtTrack; + +namespace AtTools { +class AtTrackTransformer; +} // namespace AtTools +namespace AtFITTER { +class AtFitterOld; +} // namespace AtFITTER +namespace genfit { +class Track; +} // namespace genfit + +class [[deprecated]] AtFitterTaskOld : public FairTask { + +public: + // AtFitterTask(); + ~AtFitterTaskOld() = default; + AtFitterTaskOld(std::unique_ptr fitter); + + void SetInputBranch(TString branchName); + void SetOutputBranch(TString branchName); + void SetPersistence(Bool_t value = kTRUE); + + virtual InitStatus Init(); + virtual void SetParContainers(); + virtual void Exec(Option_t *opt); + +private: + TString fInputBranchName; + TString fOutputBranchName; + + Bool_t fIsPersistence; //!< Persistence check variable + + std::unique_ptr fFitter; + AtDigiPar *fPar{nullptr}; + TClonesArray *fPatternEventArray; + TClonesArray fTrackingEventArray; + + std::size_t fEventCnt{0}; + + ClassDef(AtFitterTaskOld, 1); +}; + +#endif diff --git a/AtReconstruction/AtReconstructionLinkDef.h b/AtReconstruction/AtReconstructionLinkDef.h index 1526e0531..bd4887400 100755 --- a/AtReconstruction/AtReconstructionLinkDef.h +++ b/AtReconstruction/AtReconstructionLinkDef.h @@ -46,20 +46,25 @@ #pragma link C++ namespace kf; #pragma link C++ namespace kf::util; +#pragma link C++ namespace AtFITTER; +#pragma link C++ class AtFitterTask + ; +#pragma link C++ class AtFITTER::AtFitter - ; + /* Classes that depend on Genfit2 */ #pragma link C++ class genfit::AtSpacepointMeasurement + ; -#pragma link C++ class AtFITTER::AtFitter + ; +#pragma link C++ class AtFITTER::AtFitterOld + ; +#pragma link C++ class AtFITTER::AtFitterTaskOld + ; #pragma link C++ class AtFITTER::AtGenfit + ; -#pragma link C++ namespace AtFITTER; -#pragma link C++ class AtFitterTask + ; #pragma link C++ namespace MCFitter; #pragma link C++ class MCFitter::AtParameterDistribution - !; #pragma link C++ class MCFitter::AtUniformDistribution - !; #pragma link C++ class MCFitter::AtStudentDistribution - !; #pragma link C++ class MCFitter::AtMCFitter - !; +#pragma link C++ class MCFitter::AtMCFitterOld - !; #pragma link C++ class MCFitter::AtMCFission - !; #pragma link C++ class AtMCFitterTask + ; +#pragma link C++ class AtMCFitterTaskOld + ; /* Tasks in AtReconstruction */ #pragma link C++ class AtPSAtask + ; diff --git a/AtReconstruction/CMakeLists.txt b/AtReconstruction/CMakeLists.txt index a4bc5296b..cc40e6d9d 100755 --- a/AtReconstruction/CMakeLists.txt +++ b/AtReconstruction/CMakeLists.txt @@ -111,9 +111,16 @@ set(SRCS AtFitter/ParameterDistributions/AtUniformDistribution.cxx AtFitter/ParameterDistributions/AtStudentDistribution.cxx + AtFitter/AtFitter.cxx AtFitter/AtMCFitter.cxx AtFitter/AtMCFitterTask.cxx AtFitter/AtMCFission.cxx + AtFitterTask.cxx + + # Deprecated... + AtFitter/AtMCFitterOld.cxx + AtFitter/AtMCFitterTaskOld.cxx + ) @@ -153,10 +160,12 @@ endif() if(GENFIT2_FOUND) set(SRCS ${SRCS} - AtFitter/AtFitter.cxx AtFitter/AtGenfit.cxx AtFitter/AtSpacePointMeasurement.cxx - AtFitterTask.cxx + + # Deprecated, to be removed: + AtFitter/AtFitterOld.cxx + AtFitter/AtFitterTaskOld.cxx ) set(DEPENDENCIES ${DEPENDENCIES} GENFIT2::genfit2 diff --git a/AtTools/AtKinematics.cxx b/AtTools/AtKinematics.cxx index 2431d34f2..d924e5814 100644 --- a/AtTools/AtKinematics.cxx +++ b/AtTools/AtKinematics.cxx @@ -304,4 +304,13 @@ double EtoA(double mass) return mass / 931.5; } +std::tuple GetMomFromBrho(double mass, int Z, double brho) +{ + const Double_t M_Ener = mass * 931.49401; // In MeV + Double_t p = brho * Z * (2.99792458 / 10.0) * 1000.0; // In MeV + Double_t E = TMath::Sqrt(TMath::Power(p, 2) + TMath::Power(M_Ener, 2)) - M_Ener; // In MeV + + return std::make_tuple(p, E); +} + } // namespace AtTools::Kinematics diff --git a/AtTools/AtKinematics.h b/AtTools/AtKinematics.h index b0d343fd2..274dea9b0 100644 --- a/AtTools/AtKinematics.h +++ b/AtTools/AtKinematics.h @@ -73,6 +73,7 @@ double GetBeta(double p, double mass); double GetRelMom(double gamma, double mass); double AtoE(double Amu); double EtoA(double mass); +std::tuple GetMomFromBrho(double mass, int Z, double brho); template ROOT::Math::PxPyPzEVector Get4Vector(Vector mom, double m) diff --git a/AtTools/AtTrackTransformer.cxx b/AtTools/AtTrackTransformer.cxx index fb49ae267..b3c80d8c5 100644 --- a/AtTools/AtTrackTransformer.cxx +++ b/AtTools/AtTrackTransformer.cxx @@ -266,6 +266,26 @@ void AtTools::AtTrackTransformer::ClusterizeSmooth3D(AtTrack &track, Float_t rad } // if array size } +Bool_t AtTools::AtTrackTransformer::FindVertexTrack(AtTrack *trA, AtTrack *trB) +{ + // Determination of first hit distance. NB: Assuming both tracks have the same angle sign + Double_t vertexA = 0.0; + Double_t vertexB = 0.0; + if (trA->GetGeoTheta() * TMath::RadToDeg() < 90) { + auto iniClusterA = trA->GetHitClusterArray()->back(); + auto iniClusterB = trB->GetHitClusterArray()->back(); + vertexA = 1000.0 - iniClusterA.GetPosition().Z(); + vertexB = 1000.0 - iniClusterB.GetPosition().Z(); + } else if (trA->GetGeoTheta() * TMath::RadToDeg() > 90) { + auto iniClusterA = trA->GetHitClusterArray()->front(); + auto iniClusterB = trB->GetHitClusterArray()->front(); + vertexA = iniClusterA.GetPosition().Z(); + vertexB = iniClusterB.GetPosition().Z(); + } + + return vertexA < vertexB; +} + const std::tuple AtTools::AtTrackTransformer::GetPIDFromHits(AtTrack &track, Double_t theta) { @@ -310,3 +330,115 @@ const std::tuple AtTools::AtTrackTransformer::GetPIDFromHits return std::forward_as_tuple(dedx, eloss); } + +Bool_t AtTools::AtTrackTransformer::MergeTracks(std::vector *trackCandSource, + std::vector *trackDest, Bool_t enableSingleVertexTrack, + Double_t clusterRadius, Double_t clusterDistance) +{ + + Bool_t toMerge = kFALSE; + + Int_t addHitCnt = 0; + // Find the track closer to vertex + std::sort(trackCandSource->begin(), trackCandSource->end(), + [this](AtTrack *trA, AtTrack *trB) { return FindVertexTrack(trA, trB); }); + + // Track stitching from vertex + AtTrack *vertexTrack = *trackCandSource->begin(); + + if (enableSingleVertexTrack) { + + // Mark all tracks as merged + for (auto track : *trackCandSource) + track->SetIsMerged(kTRUE); + + trackDest->push_back(*vertexTrack); + return true; + } + + // Check if the candidate vertex track was merged + if (vertexTrack->GetIsMerged()) + return kFALSE; + else + vertexTrack->SetIsMerged(kTRUE); + + // If enabled, choose only the track closest to vertex (i.e. first one of the collection of candidates) + // TODO: Select by number of points + + for (auto it = trackCandSource->begin() + 1; it != trackCandSource->end(); ++it) { + // NB: These tracks were previously marked to merge. If merging fails they should be discarded. + AtTrack *trackToMerge = *(it); + toMerge = kFALSE; + + // Skip trackes flagged as merged + if (!trackToMerge->GetIsMerged()) { + trackToMerge->SetIsMerged(kTRUE); + } else + continue; + + Double_t endVertexZ = 0.0; + Double_t iniMergeZ = 0.0; + std::cout << " Trying to merge ... " + << "\n"; + std::cout << " Vertex track " << vertexTrack->GetTrackID() << " - Track to Merge " << trackToMerge->GetTrackID() + << "\n"; + // Check relative position between end and begin of each track using Hit Clusters + std::cout << " Vertex angle " << vertexTrack->GetGeoTheta() * TMath::RadToDeg() << "\n"; + if (vertexTrack->GetGeoTheta() * TMath::RadToDeg() < 90) { + auto endClusterVertex = vertexTrack->GetHitClusterArray()->front(); + auto iniClusterMerge = trackToMerge->GetHitClusterArray()->back(); + // Check separation and relative distance + endVertexZ = 1000.0 - endClusterVertex.GetPosition().Z(); + iniMergeZ = 1000.0 - iniClusterMerge.GetPosition().Z(); + + Double_t distance = std::sqrt((iniClusterMerge.GetPosition() - endClusterVertex.GetPosition()).Mag2()); + // std::cout << " Distance between tracks " << distance << "\n"; + // std::cout << " Ini Merge " << iniMergeZ << " - endVertexZ " << endVertexZ << "\n"; + if (((iniMergeZ + 10.0) > endVertexZ) && distance < 200) { + toMerge = kTRUE; + } + + } else if (vertexTrack->GetGeoTheta() * TMath::RadToDeg() > 90) { + auto endClusterVertex = vertexTrack->GetHitClusterArray()->back(); + auto iniClusterMerge = trackToMerge->GetHitClusterArray()->front(); + // Check separation and relative distance + endVertexZ = endClusterVertex.GetPosition().Z(); + iniMergeZ = iniClusterMerge.GetPosition().Z(); + Double_t distance = std::sqrt((iniClusterMerge.GetPosition() - endClusterVertex.GetPosition()).Mag2()); + // std::cout<<" Distance between tracks "< endVertexZ) && + distance < 100) { // NB: Distance between parts of the backward tracks is more critical + toMerge = kTRUE; + } + } + + if (toMerge) { + + std::cout << " --- Merging Succeeded! Vertex track " << vertexTrack->GetTrackID() << " - Track to Merge " + << trackToMerge->GetTrackID() << "\n"; + for (const auto &hit : trackToMerge->GetHitArray()) { + + vertexTrack->AddHit(hit->Clone()); // TODO: Look at code and see if this can be a move instead of a copy + ++addHitCnt; + } + + // Reclusterize after merging + vertexTrack->SortHitArrayTime(); + vertexTrack->ResetHitClusterArray(); + ClusterizeSmooth3D( + *vertexTrack, clusterRadius, + clusterDistance); // NB: It can be removed if we force reclusterization for any track in the mina program + + // TODO: Check if phi recalculatio is needed + + } else { + std::cout << " --- Merging Failed ! Vertex track " << vertexTrack->GetTrackID() << " - Track to Merge " + << trackToMerge->GetTrackID() << "\n"; + } + } + + trackDest->push_back(*vertexTrack); + + return toMerge; +} diff --git a/AtTools/AtTrackTransformer.h b/AtTools/AtTrackTransformer.h index b0c20db3c..acb93b78a 100644 --- a/AtTools/AtTrackTransformer.h +++ b/AtTools/AtTrackTransformer.h @@ -16,6 +16,11 @@ class AtTrackTransformer { void ClusterizeSmooth3D(AtTrack &track, Float_t radius, Float_t distance); const std::tuple GetPIDFromHits(AtTrack &track, Double_t theta); + Bool_t FindVertexTrack(AtTrack *trA, AtTrack *trB); + + Bool_t MergeTracks(std::vector *trackCandSource, std::vector *trackDest, + Bool_t enableSingleVertexTrack, Double_t clusterRadius, Double_t clusterDistance); + private: }; From 888d09164dbd0797a85fdd9406ab6335399644a4 Mon Sep 17 00:00:00 2001 From: RealAurio Date: Mon, 1 Sep 2025 16:52:46 +0200 Subject: [PATCH 2/9] Renamed and most changes regarding rebranding of AtFitMetadata (previously AtFitResult). --- AtData/AtDataLinkDef.h | 4 ++-- AtData/{AtFitResult.cxx => AtFitMetadata.cxx} | 0 AtData/{AtFitResult.h => AtFitMetadata.h} | 0 ...FitTrackResult.cxx => AtFitTrackMetadata.cxx} | 8 ++++---- .../{AtFitTrackResult.h => AtFitTrackMetadata.h} | 14 +++++++------- AtData/AtFittedTrack.h | 10 +++++----- AtData/AtMCResult.cxx | 2 +- AtData/AtMCResult.h | 4 ++-- AtData/CMakeLists.txt | 4 ++-- AtReconstruction/AtFitter/AtFitter.cxx | 2 +- AtReconstruction/AtFitter/AtFitter.h | 16 ++++++++-------- AtReconstruction/AtFitterTask.cxx | 15 ++++++++------- AtReconstruction/AtFitterTask.h | 8 ++++---- 13 files changed, 44 insertions(+), 43 deletions(-) rename AtData/{AtFitResult.cxx => AtFitMetadata.cxx} (100%) rename AtData/{AtFitResult.h => AtFitMetadata.h} (100%) rename AtData/{AtFitTrackResult.cxx => AtFitTrackMetadata.cxx} (62%) rename AtData/{AtFitTrackResult.h => AtFitTrackMetadata.h} (77%) diff --git a/AtData/AtDataLinkDef.h b/AtData/AtDataLinkDef.h index 013860b29..411d87117 100644 --- a/AtData/AtDataLinkDef.h +++ b/AtData/AtDataLinkDef.h @@ -43,8 +43,8 @@ #pragma link C++ enum AtPatterns::PatternType; #pragma link C++ function AtPatterns::CreatePattern; -#pragma link C++ class AtFitResult + ; -#pragma link C++ class AtFitTrackResult + ; +#pragma link C++ class AtFitMetadata + ; +#pragma link C++ class AtFitTrackMetadata + ; #pragma link C++ class MCFitter::AtMCResult + ; #pragma link C++ class AtTrackingEventOld + ; diff --git a/AtData/AtFitResult.cxx b/AtData/AtFitMetadata.cxx similarity index 100% rename from AtData/AtFitResult.cxx rename to AtData/AtFitMetadata.cxx diff --git a/AtData/AtFitResult.h b/AtData/AtFitMetadata.h similarity index 100% rename from AtData/AtFitResult.h rename to AtData/AtFitMetadata.h diff --git a/AtData/AtFitTrackResult.cxx b/AtData/AtFitTrackMetadata.cxx similarity index 62% rename from AtData/AtFitTrackResult.cxx rename to AtData/AtFitTrackMetadata.cxx index 2c93509b9..d2627dfcb 100644 --- a/AtData/AtFitTrackResult.cxx +++ b/AtData/AtFitTrackMetadata.cxx @@ -1,12 +1,12 @@ -#include "AtFitTrackResult.h" +#include "AtFitTrackMetadata.h" #include -ClassImp(AtFitTrackResult); +ClassImp(AtFitTrackMetadata); -void AtFitTrackResult::Print() const +void AtFitTrackMetadata::Print() const { - std::cout << " Fit result for track with ID " << fTrackID << ":" << std::endl; + std::cout << " Fit metadata for track with ID " << fTrackID << ":" << std::endl; std::cout << " Statistics: " << std::endl; std::cout << " PValue = " << fPValue << std::endl; diff --git a/AtData/AtFitTrackResult.h b/AtData/AtFitTrackMetadata.h similarity index 77% rename from AtData/AtFitTrackResult.h rename to AtData/AtFitTrackMetadata.h index 2960a84c9..b21ad6729 100644 --- a/AtData/AtFitTrackResult.h +++ b/AtData/AtFitTrackMetadata.h @@ -1,5 +1,5 @@ -#ifndef ATFITTRACKRESULT_H -#define ATFITTRACKRESULT_H +#ifndef ATFITTRACKMETADATA_H +#define ATFITTRACKMETADATA_H #include // for Double_t, THashConsistencyHolder, ClassDefOverride #include @@ -11,9 +11,9 @@ class TMemberInspector; /** * Class for storing the result of the fit of an AtTrack from an AtFitter class. */ -class AtFitTrackResult : public TObject { +class AtFitTrackMetadata : public TObject { protected: - // Copy of the statistics parameters which are also saved in AtFittedTrack. + // Statistics parameters of the fit. Double_t fPValue{0}; Double_t fChi2{0}; Int_t fNdf{0}; @@ -23,8 +23,8 @@ class AtFitTrackResult : public TObject { Int_t fTrackID{-1}; public: - AtFitTrackResult() = default; - ~AtFitTrackResult() = default; + AtFitTrackMetadata() = default; + ~AtFitTrackMetadata() = default; void SetPValue(Double_t value) { fPValue = value; } void SetChi2(Double_t value) { fChi2 = value; } @@ -40,7 +40,7 @@ class AtFitTrackResult : public TObject { virtual void Print() const; - ClassDefOverride(AtFitTrackResult, 1); + ClassDefOverride(AtFitTrackMetadata, 1); }; #endif diff --git a/AtData/AtFittedTrack.h b/AtData/AtFittedTrack.h index fce53b138..b5bb0d4d8 100644 --- a/AtData/AtFittedTrack.h +++ b/AtData/AtFittedTrack.h @@ -1,7 +1,7 @@ #ifndef ATFITTEDTRACK_H #define ATFITTEDTRACK_H -#include "AtFitTrackResult.h" +#include "AtFitTrackMetadata.h" #include #include @@ -25,7 +25,7 @@ class TMemberInspector; class AtFittedTrack : public TObject { public: using XYZVector = ROOT::Math::XYZVector; - using TrackResultPtr = std::unique_ptr; + using TrackMetadataPtr = std::unique_ptr; struct Kinematics { Double_t kineticEnergy{-1}; // Kinetic energy @@ -81,7 +81,7 @@ class AtFittedTrack : public TObject { Statistics fStats; // Copy of the AtFitTrackResult object corresponding to the fit used for this track. - TrackResultPtr fTrackResult{nullptr}; + TrackMetadataPtr fTrackMetadata{nullptr}; public: AtFittedTrack() = default; @@ -128,7 +128,7 @@ class AtFittedTrack : public TObject { fStats.fitConverged = conv; } - void SetTrackResult(TrackResultPtr trackResult) { fTrackResult = std::move(trackResult); } + void SetTrackMetadata(TrackMetadataPtr trackMetadata) { fTrackMetadata = std::move(trackMetadata); } const Int_t GetTrackID() { return fTrackID; } @@ -143,7 +143,7 @@ class AtFittedTrack : public TObject { const TrackProperties GetTrackProperties() { return fTrackProperties; } const Statistics GetStats() { return fStats; } - TrackResultPtr &GetTrackResult() { return fTrackResult; } + TrackMetadataPtr &GetTrackMetadata() { return fTrackMetadata; } ClassDef(AtFittedTrack, 2); }; diff --git a/AtData/AtMCResult.cxx b/AtData/AtMCResult.cxx index 053e926fc..23de778d3 100644 --- a/AtData/AtMCResult.cxx +++ b/AtData/AtMCResult.cxx @@ -7,7 +7,7 @@ namespace MCFitter { void AtMCResult::Print() const { - AtFitTrackResult::Print(); + AtFitTrackMetadata::Print(); std::cout << " MC fit specifics: " << std::endl; std::cout << " Iteration = " << fIterNum << std::endl; diff --git a/AtData/AtMCResult.h b/AtData/AtMCResult.h index d4d48246b..98c28c405 100644 --- a/AtData/AtMCResult.h +++ b/AtData/AtMCResult.h @@ -1,7 +1,7 @@ #ifndef ATMCRESULT_H #define ATMCRESULT_H -#include "AtFitTrackResult.h" +#include "AtFitTrackMetadata.h" #include // for Double_t, THashConsistencyHolder, ClassDefOverride #include @@ -17,7 +17,7 @@ namespace MCFitter { /** * Class for storing the result of an iteration in the AtMCFitter method. */ -class AtMCResult : public AtFitTrackResult { +class AtMCResult : public AtFitTrackMetadata { public: using ParamMap = std::map; diff --git a/AtData/CMakeLists.txt b/AtData/CMakeLists.txt index b3c0a5d1a..3ab791cff 100644 --- a/AtData/CMakeLists.txt +++ b/AtData/CMakeLists.txt @@ -54,8 +54,8 @@ set(SRCS AtPattern/AtPadPlaneElement.cxx AtFittedTrack.cxx - AtFitResult.cxx - AtFitTrackResult.cxx + AtFitMetadata.cxx + AtFitTrackMetadata.cxx AtMCResult.cxx diff --git a/AtReconstruction/AtFitter/AtFitter.cxx b/AtReconstruction/AtFitter/AtFitter.cxx index 9cca3032d..fc6d7ade1 100644 --- a/AtReconstruction/AtFitter/AtFitter.cxx +++ b/AtReconstruction/AtFitter/AtFitter.cxx @@ -5,7 +5,7 @@ ClassImp(AtFITTER::AtFitter); void AtFITTER::AtFitter::Reset() { fPatternEvent = nullptr; - fFitResult = nullptr; + fFitMetadata = nullptr; fRawEvent = nullptr; fEvent = nullptr; } diff --git a/AtReconstruction/AtFitter/AtFitter.h b/AtReconstruction/AtFitter/AtFitter.h index 32233e28b..a745550eb 100644 --- a/AtReconstruction/AtFitter/AtFitter.h +++ b/AtReconstruction/AtFitter/AtFitter.h @@ -1,7 +1,7 @@ #ifndef AtFITTER_H #define AtFITTER_H -#include "AtFitResult.h" +#include "AtFitMetadata.h" #include #include @@ -26,17 +26,17 @@ namespace AtFITTER { class AtFitter : public TObject { public: - using TrackResultPtr = std::unique_ptr; - using TrackResultsVector = std::vector; - using TrackResultsSet = - std::set>; + using TrackMetadataPtr = std::unique_ptr; + using TrackMetadatasVector = std::vector; + using TrackMetadatasSet = + std::set>; protected: // Pointer to the AtPatternEvent to be fitted. AtPatternEvent *fPatternEvent{nullptr}; // Pointer to the AtFitResult object in which store the fit metadata. - AtFitResult *fFitResult{nullptr}; + AtFitMetadata *fFitMetadata{nullptr}; // Pointers to AtRawEvent and AtEvent. In case some specific fitter needs to access information // in any of those branches. @@ -44,7 +44,7 @@ class AtFitter : public TObject { AtEvent *fEvent{nullptr}; // Compare function that will be used to sort the fit results for a given track. - virtual bool CompareTrackFitsFunction(const TrackResultPtr &trackResultA, const TrackResultPtr &trackResultB) = 0; + virtual bool CompareTrackFitsFunction(const TrackMetadataPtr &trackMetadataA, const TrackMetadataPtr &trackMetadataB) = 0; public: AtFitter() = default; @@ -55,7 +55,7 @@ class AtFitter : public TObject { // Mandatory to set. void SetPatternEvent(AtPatternEvent *patternEvent) { fPatternEvent = patternEvent; } - void SetFitResult(AtFitResult *fitResult) { fFitResult = fitResult; } + void SetFitMetadata(AtFitMetadata *fitMetadata) { fFitMetadata = fitMetadata; } // Optional to set. void SetRawEvent(AtRawEvent *rawEvent) { fRawEvent = rawEvent; } diff --git a/AtReconstruction/AtFitterTask.cxx b/AtReconstruction/AtFitterTask.cxx index 5f73e967f..aaf52d044 100644 --- a/AtReconstruction/AtFitterTask.cxx +++ b/AtReconstruction/AtFitterTask.cxx @@ -2,6 +2,7 @@ #include "AtDigiPar.h" #include "AtEvent.h" +#include "AtFitMetadata.h" #include "AtFitter.h" #include "AtParsers.h" #include "AtPatternEvent.h" @@ -28,7 +29,7 @@ ClassImp(AtFitterTask); AtFitterTask::AtFitterTask(std::unique_ptr fitter) : fInputBranchName("AtPatternEvent"), fOutputBranchName("AtTrackingEvent"), fIsPersistence(kFALSE), fTrackingEventArray(TClonesArray("AtTrackingEvent", 1)), fFitter(std::move(fitter)), fRawEventBranchName(""), - fEventBranchName(""), fFitResultBranchName("") + fEventBranchName(""), fFitMetadataBranchName("") { } @@ -57,10 +58,10 @@ void AtFitterTask::SetEventBranch(TString branchName) fEventBranchName = branchName; } -void AtFitterTask::SetFitResultBranch(TString branchName) +void AtFitterTask::SetFitMetadataBranch(TString branchName) { - fFitResultBranchName = branchName; - fSaveFitResult = true; + fFitMetadataBranchName = branchName; + fSaveFitMetadata = true; } InitStatus AtFitterTask::Init() @@ -78,7 +79,7 @@ InitStatus AtFitterTask::Init() } ioMan->Register(fOutputBranchName, "AtTPC", &fTrackingEventArray, fIsPersistence); - ioMan->Register(fFitResultBranchName, "AtTPC", &fFitResultArray, fIsPersistence && fSaveFitResult); + ioMan->Register(fFitMetadataBranchName, "AtTPC", &fFitMetadataArray, fIsPersistence && fSaveFitMetadata); fRawEventArray = dynamic_cast(ioMan->GetObject(fRawEventBranchName)); if (fRawEventArray == nullptr) { @@ -132,7 +133,7 @@ void AtFitterTask::Exec(Option_t *option) fTrackingEventArray.Delete(); auto trackingEvent = dynamic_cast(fTrackingEventArray.ConstructedAt(0)); - auto fitResult = dynamic_cast(fFitResultArray.ConstructedAt(0)); + auto fitMetadata = dynamic_cast(fFitMetadataArray.ConstructedAt(0)); LOG(info) << " AtFitterTask::Exec() : Fitting event " << fEventCnt; @@ -141,7 +142,7 @@ void AtFitterTask::Exec(Option_t *option) LOG(info) << " AtFitterTask::Exec() : Number of candidate tracks : " << tracks.size(); fFitter->SetPatternEvent(patternEvent); - fFitter->SetFitResult(fitResult); + fFitter->SetFitMetadata(fitMetadata); auto fittedTracks = fFitter->ProcessEvent(); LOG(info) << " AtFitterTask::Exec() : Number of fitted tracks : " << fittedTracks.size(); diff --git a/AtReconstruction/AtFitterTask.h b/AtReconstruction/AtFitterTask.h index 74a706227..c950382c0 100644 --- a/AtReconstruction/AtFitterTask.h +++ b/AtReconstruction/AtFitterTask.h @@ -54,7 +54,7 @@ class AtFitterTask : public FairTask { void SetRawEventBranch(TString branchName); void SetEventBranch(TString branchName); - void SetFitResultBranch(TString branchName); + void SetFitMetadataBranch(TString branchName); virtual InitStatus Init(); virtual void SetParContainers(); @@ -80,9 +80,9 @@ class AtFitterTask : public FairTask { TClonesArray *fEventArray; // Include the option to store all the fit metadata in a AtFitResult branch. - TString fFitResultBranchName; - TClonesArray fFitResultArray; - bool fSaveFitResult{false}; + TString fFitMetadataBranchName; + TClonesArray fFitMetadataArray; + bool fSaveFitMetadata{false}; ClassDef(AtFitterTask, 2); }; From 43f4ae3050d6c88ec3f3f4b8eba3d708293e5ecb Mon Sep 17 00:00:00 2001 From: RealAurio Date: Mon, 1 Sep 2025 16:53:20 +0200 Subject: [PATCH 3/9] Last changes, in different commit due to renaming of the file. --- AtData/AtFitMetadata.cxx | 4 ++-- AtData/AtFitMetadata.h | 26 +++++++++++++------------- 2 files changed, 15 insertions(+), 15 deletions(-) diff --git a/AtData/AtFitMetadata.cxx b/AtData/AtFitMetadata.cxx index 5412b7864..8e6906a91 100644 --- a/AtData/AtFitMetadata.cxx +++ b/AtData/AtFitMetadata.cxx @@ -1,3 +1,3 @@ -#include "AtFitResult.h" +#include "AtFitMetadata.h" -ClassImp(AtFitResult); +ClassImp(AtFitMetadata); diff --git a/AtData/AtFitMetadata.h b/AtData/AtFitMetadata.h index 5542d29ff..597f10da9 100644 --- a/AtData/AtFitMetadata.h +++ b/AtData/AtFitMetadata.h @@ -1,7 +1,7 @@ -#ifndef ATFITRESULT_H -#define ATFITRESULT_H +#ifndef ATFITMETADATA_H +#define ATFITMETADATA_H -#include "AtFitTrackResult.h" +#include "AtFitTrackMetadata.h" #include @@ -20,30 +20,30 @@ class TMemberInspector; /** * Class for storing the result of the fit for the entire AtTrackingEvent from an AtFitter class. */ -class AtFitResult : public TObject { +class AtFitMetadata : public TObject { public: - using TrackResultPtr = std::unique_ptr; - using TrackResultsVector = std::vector; - using ResultsMap = std::map; + using TrackMetadataPtr = std::unique_ptr; + using TrackMetadatasVector = std::vector; + using MetadatasMap = std::map; protected: // Vector to store the results for all different fits done to all tracks in the event. - ResultsMap fResults; + MetadatasMap fMetadatas; // Event ID for which this fit was done. ULong_t fEventID; public: - AtFitResult() = default; - ~AtFitResult() = default; + AtFitMetadata() = default; + ~AtFitMetadata() = default; - void SetTrackResultsVector(Int_t trackID, TrackResultsVector results) { fResults[trackID] = std::move(results); } + void SetTrackMetadatasVector(Int_t trackID, TrackMetadatasVector metadatas) { fMetadatas[trackID] = std::move(metadatas); } void SetEventID(ULong_t id) { fEventID = id; } - TrackResultsVector &GetTrackResultsVector(Int_t trackID) { return fResults.at(trackID); } + TrackMetadatasVector &GetTrackMetadatasVector(Int_t trackID) { return fMetadatas.at(trackID); } - ClassDefOverride(AtFitResult, 1); + ClassDefOverride(AtFitMetadata, 1); }; #endif From d89796fc2839ec2ce5ce4b37a8db35e4edb41d00 Mon Sep 17 00:00:00 2001 From: RealAurio Date: Mon, 1 Sep 2025 17:23:03 +0200 Subject: [PATCH 4/9] Changes some comments, removed statistics from AtFittedTrack and clang-formatted. --- AtData/AtFitMetadata.h | 13 +++++++++++-- AtData/AtFittedTrack.h | 25 ++----------------------- AtReconstruction/AtFitter/AtFitter.h | 3 ++- 3 files changed, 15 insertions(+), 26 deletions(-) diff --git a/AtData/AtFitMetadata.h b/AtData/AtFitMetadata.h index 597f10da9..9109e0e38 100644 --- a/AtData/AtFitMetadata.h +++ b/AtData/AtFitMetadata.h @@ -27,7 +27,13 @@ class AtFitMetadata : public TObject { using MetadatasMap = std::map; protected: - // Vector to store the results for all different fits done to all tracks in the event. + /** + * Map to store the metadatas for all different fits done to all tracks in the event. + * The Int_t corresponds to the trackID for which the metadatas correspond to. + * The vector of AtFitTrackMetadata contains the different metadatas for all fits + * that have been done for the track (for example, different assumptions for the + * particles of the track, different initial conditions, etc...). + */ MetadatasMap fMetadatas; // Event ID for which this fit was done. @@ -37,7 +43,10 @@ class AtFitMetadata : public TObject { AtFitMetadata() = default; ~AtFitMetadata() = default; - void SetTrackMetadatasVector(Int_t trackID, TrackMetadatasVector metadatas) { fMetadatas[trackID] = std::move(metadatas); } + void SetTrackMetadatasVector(Int_t trackID, TrackMetadatasVector metadatas) + { + fMetadatas[trackID] = std::move(metadatas); + } void SetEventID(ULong_t id) { fEventID = id; } diff --git a/AtData/AtFittedTrack.h b/AtData/AtFittedTrack.h index b5bb0d4d8..407cfa3c9 100644 --- a/AtData/AtFittedTrack.h +++ b/AtData/AtFittedTrack.h @@ -36,18 +36,10 @@ class AtFittedTrack : public TObject { Double_t phiXtr{-1}; // Extrapolated phi scattering angle }; - struct Statistics { - Double_t pValue{-1}; // Probability for rejecting the fit hypothesis - Double_t chi2{-1}; // Chi2 of the fit - Int_t NDF{-1}; // Number of degrees of freedom for the fit - Double_t redChi2{-1}; // Reduced chi2 - Bool_t fitConverged{0}; // Whether or not the fit managed to converge - }; - struct ParticleInfo { TString idPDG{""}; // PDG code of the particle - Int_t charge{0}; // Charge of the particle - Double_t mass{-1}; // Mass of the particle + Int_t charge{0}; // Charge number of the particle + Double_t mass{-1}; // Mass of the particle in amu }; struct TrackProperties { @@ -77,9 +69,6 @@ class AtFittedTrack : public TObject { // Track properties. TrackProperties fTrackProperties; - // Parameters regarding the statistics of the fit. - Statistics fStats; - // Copy of the AtFitTrackResult object corresponding to the fit used for this track. TrackMetadataPtr fTrackMetadata{nullptr}; @@ -119,15 +108,6 @@ class AtFittedTrack : public TObject { fTrackProperties.trackPoints = trackPoints; } - void SetStats(Double_t pvalue, Double_t chi2, Int_t ndf, Bool_t conv) - { - fStats.pValue = pvalue; - fStats.chi2 = chi2; - fStats.NDF = ndf; - fStats.redChi2 = chi2 / ndf; - fStats.fitConverged = conv; - } - void SetTrackMetadata(TrackMetadataPtr trackMetadata) { fTrackMetadata = std::move(trackMetadata); } const Int_t GetTrackID() { return fTrackID; } @@ -141,7 +121,6 @@ class AtFittedTrack : public TObject { const XYZVector GetVertex() { return fVertex[0]; } const TrackProperties GetTrackProperties() { return fTrackProperties; } - const Statistics GetStats() { return fStats; } TrackMetadataPtr &GetTrackMetadata() { return fTrackMetadata; } diff --git a/AtReconstruction/AtFitter/AtFitter.h b/AtReconstruction/AtFitter/AtFitter.h index a745550eb..4414388fc 100644 --- a/AtReconstruction/AtFitter/AtFitter.h +++ b/AtReconstruction/AtFitter/AtFitter.h @@ -44,7 +44,8 @@ class AtFitter : public TObject { AtEvent *fEvent{nullptr}; // Compare function that will be used to sort the fit results for a given track. - virtual bool CompareTrackFitsFunction(const TrackMetadataPtr &trackMetadataA, const TrackMetadataPtr &trackMetadataB) = 0; + virtual bool + CompareTrackFitsFunction(const TrackMetadataPtr &trackMetadataA, const TrackMetadataPtr &trackMetadataB) = 0; public: AtFitter() = default; From ef876d7172a7d74c5a052bd8a08bbbaa1f5685fa Mon Sep 17 00:00:00 2001 From: RealAurio Date: Tue, 2 Sep 2025 17:51:48 +0200 Subject: [PATCH 5/9] Brought back the old getters and setters, as well as split the Kinematics in 2. --- AtData/AtFittedTrack.cxx | 18 +++-- AtData/AtFittedTrack.h | 169 ++++++++++++++++++++++++++++++++++----- 2 files changed, 160 insertions(+), 27 deletions(-) diff --git a/AtData/AtFittedTrack.cxx b/AtData/AtFittedTrack.cxx index 9a5756bc2..940a59002 100644 --- a/AtData/AtFittedTrack.cxx +++ b/AtData/AtFittedTrack.cxx @@ -2,8 +2,7 @@ ClassImp(AtFittedTrack); -void AtFittedTrack::SetKinematics(int particleIdx, Double_t energy, Double_t theta, Double_t phi, Double_t energyxtr, - Double_t thetaxtr, Double_t phixtr) +void AtFittedTrack::SetKinematics(int particleIdx, Double_t energy, Double_t theta, Double_t phi) { while (particleIdx >= fKinematics.size()) { Kinematics newKinematics; @@ -13,9 +12,18 @@ void AtFittedTrack::SetKinematics(int particleIdx, Double_t energy, Double_t the fKinematics[particleIdx].kineticEnergy = energy; fKinematics[particleIdx].theta = theta; fKinematics[particleIdx].phi = phi; - fKinematics[particleIdx].kineticEnergyXtr = energyxtr; - fKinematics[particleIdx].thetaXtr = thetaxtr; - fKinematics[particleIdx].phiXtr = phixtr; +} + +void AtFittedTrack::SetKinematicsXtr(int particleIdx, Double_t energyxtr, Double_t thetaxtr, Double_t phixtr) +{ + while (particleIdx >= fKinematicsXtr.size()) { + Kinematics newKinematics; + fKinematicsXtr.push_back(newKinematics); + } + + fKinematicsXtr[particleIdx].kineticEnergy = energyxtr; + fKinematicsXtr[particleIdx].theta = thetaxtr; + fKinematicsXtr[particleIdx].phi = phixtr; } void AtFittedTrack::SetParticleInfo(int particleIdx, std::string pdg, Int_t charge, Double_t mass) diff --git a/AtData/AtFittedTrack.h b/AtData/AtFittedTrack.h index 407cfa3c9..62727a81d 100644 --- a/AtData/AtFittedTrack.h +++ b/AtData/AtFittedTrack.h @@ -28,12 +28,9 @@ class AtFittedTrack : public TObject { using TrackMetadataPtr = std::unique_ptr; struct Kinematics { - Double_t kineticEnergy{-1}; // Kinetic energy - Double_t theta{-1}; // Theta scattering angle - Double_t phi{-1}; // Phi scattering angle - Double_t kineticEnergyXtr{-1}; // Extrapolated kinetic energy - Double_t thetaXtr{-1}; // Extrapolated theta scattering angle - Double_t phiXtr{-1}; // Extrapolated phi scattering angle + Double_t kineticEnergy{-1}; // Kinetic energy + Double_t theta{-1}; // Theta scattering angle + Double_t phi{-1}; // Phi scattering angle }; struct ParticleInfo { @@ -59,6 +56,7 @@ class AtFittedTrack : public TObject { // Kinematic variables obtained by the fit. std::vector fKinematics; + std::vector fKinematicsXtr; // Particle information. std::vector fParticleInfo; @@ -72,30 +70,45 @@ class AtFittedTrack : public TObject { // Copy of the AtFitTrackResult object corresponding to the fit used for this track. TrackMetadataPtr fTrackMetadata{nullptr}; + // Deprecated members needed to keep support of deprecated methods. + [[deprecated("No PRA information is supposed to be kept in AtFittedTrack.")]] Double_t fEnergyPRA{0}; + [[deprecated("No PRA information is supposed to be kept in AtFittedTrack.")]] Double_t fThetaPRA{0}; + [[deprecated("No PRA information is supposed to be kept in AtFittedTrack.")]] Double_t fPhiPRA{0}; + [[deprecated("No PRA information is supposed to be kept in AtFittedTrack.")]] XYZVector fInitialPosPRA; + [[deprecated("No backwards statistics is supposed to be stored. May not be needed in all fitters.")]] Float_t fBChi2{ + 0}; + [[deprecated("No backwards statistics is supposed to be stored. May not be needed in all fitters.")]] Float_t fBNdf{ + 0}; + [[deprecated("Not all experiments are done with magnetic field.")]] Float_t fBrho{0}; + [[deprecated("No IC information is supposed to be kept in AtFittedTrack.")]] Float_t fIonChamberEnergy{0}; + [[deprecated("No IC information is supposed to be kept in AtFittedTrack.")]] Int_t fIonChamberTime{0}; + [[deprecated("No excitation energy is supposed to be kept in AtFittedTrack.")]] Float_t fExcitationEnergy{0}; + [[deprecated("No excitation energy is supposed to be kept in AtFittedTrack.")]] Float_t fExcitationEnergyXtr{0}; + public: AtFittedTrack() = default; ~AtFittedTrack() = default; void SetTrackID(Int_t trackid) { fTrackID = trackid; } - void SetKinematics(int particleIdx, Double_t energy, Double_t theta, Double_t phi, Double_t energyxtr, - Double_t thetaxtr, Double_t phixtr); + void SetKinematics(int particleIdx, Double_t energy, Double_t theta, Double_t phi); + void SetKinematicsXtr(int particleIdx, Double_t energyxtr, Double_t thetaxtr, Double_t phixtr); void SetParticleInfo(int particleIdx, std::string pdg, Int_t charge, Double_t mass); void SetVertex(int particleIdx, XYZVector point); - void - SetKinematics(Double_t energy, Double_t theta, Double_t phi, Double_t energyxtr, Double_t thetaxtr, Double_t phixtr) + void SetKinematics(Double_t energy, Double_t theta, Double_t phi) { SetKinematics(0, energy, theta, phi); } + void SetKinematicsXtr(Double_t energyxtr, Double_t thetaxtr, Double_t phixtr) { - SetKinematics(0, energy, theta, phi, energyxtr, thetaxtr, phixtr); + SetKinematicsXtr(0, energyxtr, thetaxtr, phixtr); } void SetParticleInfo(std::string pdg, Int_t charge, Double_t mass) { SetParticleInfo(0, pdg, charge, mass); } void SetVertex(XYZVector point) { SetVertex(0, point); } - void SetTrackProperties(XYZVector initialPosition, XYZVector initialPositionXtr, Double_t extrapolatedDistance, - Double_t distancePOCA, Double_t trackLength, Double_t trackLengthXtr, - Double_t estimateTotalCharge, Int_t trackPoints) + void SetTrackPropertiesStruct(XYZVector initialPosition, XYZVector initialPositionXtr, Double_t extrapolatedDistance, + Double_t distancePOCA, Double_t trackLength, Double_t trackLengthXtr, + Double_t estimateTotalCharge, Int_t trackPoints) { fTrackProperties.initialPosition = initialPosition; fTrackProperties.initialPositionXtr = initialPositionXtr; @@ -112,18 +125,130 @@ class AtFittedTrack : public TObject { const Int_t GetTrackID() { return fTrackID; } - const Kinematics GetKinematics(int particleIdx) { return fKinematics[particleIdx]; } - const ParticleInfo GetParticleInfo(int particleIdx) { return fParticleInfo[particleIdx]; } - const XYZVector GetVertex(int particleIdx) { return fVertex[particleIdx]; } - - const Kinematics GetKinematics() { return fKinematics[0]; } - const ParticleInfo GetParticleInfo() { return fParticleInfo[0]; } - const XYZVector GetVertex() { return fVertex[0]; } + const Kinematics GetKinematics(int particleIdx = 0) { return fKinematics[particleIdx]; } + const Kinematics GetKinematicsXtr(int particleIdx = 0) { return fKinematicsXtr[particleIdx]; } + const ParticleInfo GetParticleInfo(int particleIdx = 0) { return fParticleInfo[particleIdx]; } + const XYZVector GetVertex(int particleIdx = 0) { return fVertex[particleIdx]; } - const TrackProperties GetTrackProperties() { return fTrackProperties; } + const TrackProperties GetTrackPropertiesStruct() { return fTrackProperties; } TrackMetadataPtr &GetTrackMetadata() { return fTrackMetadata; } + // Old deprecated methods. + [[deprecated("Replaced by SetKinematics() and SetKinematicsXtr().")]] void + SetEnergyAngles(Float_t energy, Float_t energyxtr, Float_t theta, Float_t phi, Float_t energypra, Float_t thetapra, + Float_t phipra) + { + fKinematics[0].kineticEnergy = energy; + fKinematics[0].theta = theta; + fKinematics[0].phi = phi; + fKinematicsXtr[0].kineticEnergy = energyxtr; + fEnergyPRA = energypra; + fThetaPRA = thetapra; + fPhiPRA = phipra; + } + [[deprecated("This information now lives in the TrackProperties struct. Check the new SetTrackProperties().")]] void + SetVertexPosition(XYZVector inipos, XYZVector iniposPRA, XYZVector iniposxtr) + { + fTrackProperties.initialPosition = inipos; + fInitialPosPRA = iniposPRA; + fVertex[0] = iniposxtr; + } + [[deprecated("Statistics now live inside the AtFitTrackMetadata. Check SetTrackMetadata().")]] void + SetStats(Float_t pvalue, Float_t chi2, Float_t bchi2, Float_t ndf, Float_t bndf, Bool_t conv) + { + if (fTrackMetadata == nullptr) + fTrackMetadata = std::make_unique(); + fTrackMetadata->SetPValue(pvalue); + fTrackMetadata->SetChi2(chi2); + fTrackMetadata->SetNdf(ndf); + fTrackMetadata->SetFitConverged(conv); + fTrackMetadata->SetTrackID(fTrackID); + + fBChi2 = bchi2; + fBNdf = bndf; + } + [[deprecated("The TrackProperties have changed. Please check the new SetTrackProperties() method.")]] void + SetTrackProperties(Int_t charge, Float_t brho, Float_t eloss, Float_t dedx, std::string pdg, Int_t points) + { + fParticleInfo[0].charge = charge; + fBrho = brho; + fTrackProperties.estimateTotalCharge = eloss; + fTrackProperties.estimateDeDx = dedx; + fParticleInfo[0].idPDG = pdg; + fTrackProperties.trackPoints = points; + } + [[deprecated("Ion chamber information is no longer saved in the AtFittedTrack. This has been saved here still, but " + "refrain from using this further.")]] void + SetIonChamber(Float_t icenergy, Int_t ictime) + { + fIonChamberEnergy = icenergy; + fIonChamberTime = ictime; + } + [[deprecated("Excitation energy is no longer saved in the AtFittedTrack. This has been saved here still, but " + "refrain from using this further.")]] void + SetExcitationEnergy(Float_t exenergy, Float_t exenergyxtr) + { + fExcitationEnergy = exenergy; + fExcitationEnergyXtr = exenergyxtr; + } + [[deprecated("This information now lives in the TrackProperties struct. Check the new SetTrackProperties().")]] void + SetDistances(Float_t distancextr, Float_t length, Float_t poca) + { + fTrackProperties.extrapolatedDistance = distancextr; + fTrackProperties.trackLength = length; + fTrackProperties.distancePOCA = poca; + } + + [[deprecated("Please check GetKinematics() and GetKinematicsXtr().")]] const std::tuple< + Float_t, Float_t, Float_t, Float_t, Float_t, Float_t, Float_t> + GetEnergyAngles() + { + return std::forward_as_tuple(fKinematics[0].kineticEnergy, fKinematicsXtr[0].kineticEnergy, fKinematics[0].theta, + fKinematics[0].phi, fEnergyPRA, fThetaPRA, fPhiPRA); + } + [[deprecated("Please check GetVertex().")]] const std::tuple GetVertices() + { + return std::forward_as_tuple(fTrackProperties.initialPosition, fInitialPosPRA, fVertex[0]); + } + [[deprecated("Statistics now live inside the AtFitTrackMetadata. Check GetTrackMetadata().")]] const std::tuple< + Float_t, Float_t, Float_t, Float_t, Float_t, Bool_t> + GetStats() + { + if (fTrackMetadata == nullptr) + return std::forward_as_tuple(0, 0, 0, 0, 0, 0); + return std::forward_as_tuple(fTrackMetadata->GetPValue(), fTrackMetadata->GetChi2(), fBChi2, + fTrackMetadata->GetNdf(), fBNdf, fTrackMetadata->GetFitConverged()); + } + [[deprecated("The TrackProperties have changed. Please check the new GetTrackProperties() method.")]] const std:: + tuple + GetTrackProperties() + { + return std::forward_as_tuple(fParticleInfo[0].charge, fBrho, fTrackProperties.estimateTotalCharge, + fTrackProperties.estimateDeDx, fParticleInfo[0].idPDG.Data(), + fTrackProperties.trackPoints); + } + [[deprecated("Ion chamber information is no longer saved in the AtFittedTrack. This has been saved here still, but " + "refrain from using this further.")]] const std::tuple + GetIonChamber() + { + return std::forward_as_tuple(fIonChamberEnergy, fIonChamberTime); + } + [[deprecated("Excitation energy is no longer saved in the AtFittedTrack. This has been saved here still, but " + "refrain from using this further.")]] const std::tuple + GetExcitationEnergy() + { + return std::forward_as_tuple(fExcitationEnergy, fExcitationEnergyXtr); + } + [[deprecated( + "This information now lives in the TrackProperties struct. Check the new SetTrackProperties().")]] const std:: + tuple + GetDistances() + { + return std::forward_as_tuple(fTrackProperties.extrapolatedDistance, fTrackProperties.trackLength, + fTrackProperties.distancePOCA); + } + ClassDef(AtFittedTrack, 2); }; From e4b9e8914e32f18d461eff95b752255d7506207c Mon Sep 17 00:00:00 2001 From: RealAurio Date: Wed, 3 Sep 2025 12:25:51 +0200 Subject: [PATCH 6/9] Changed the interface of AtFitter, changed the name of namespace from AtFITTER to EventFit. --- AtReconstruction/AtFitter/AtFitter.cxx | 23 +++++++++--- AtReconstruction/AtFitter/AtFitter.h | 43 +++++++--------------- AtReconstruction/AtFitter/AtGenfit.h | 6 ++- AtReconstruction/AtFitterTask.cxx | 35 +++++++----------- AtReconstruction/AtFitterTask.h | 8 ++-- AtReconstruction/AtReconstructionLinkDef.h | 7 ++-- 6 files changed, 57 insertions(+), 65 deletions(-) diff --git a/AtReconstruction/AtFitter/AtFitter.cxx b/AtReconstruction/AtFitter/AtFitter.cxx index fc6d7ade1..81e02ef4d 100644 --- a/AtReconstruction/AtFitter/AtFitter.cxx +++ b/AtReconstruction/AtFitter/AtFitter.cxx @@ -1,11 +1,22 @@ #include "AtFitter.h" -ClassImp(AtFITTER::AtFitter); +#include "AtPatternEvent.h" +#include "AtTrackingEvent.h" -void AtFITTER::AtFitter::Reset() +ClassImp(EventFit::AtFitter); + +void EventFit::AtFitter::FitEvent(AtTrackingEvent *trackingEvent, AtPatternEvent *patternEvent, + AtFitMetadata *fitMetadata, AtRawEvent *rawEvent, AtEvent *event) { - fPatternEvent = nullptr; - fFitMetadata = nullptr; - fRawEvent = nullptr; - fEvent = nullptr; + // Extract the candidate AtTracks. + std::vector tracks = patternEvent->GetTrackCand(); + + // Save the original AtTracks to the AtTrackingEvent. + trackingEvent->SetTrackArray(&tracks); + + // Iterate over the AtTracks and store the AtFittedTracks in the AtTrackingEvent. + for (auto track : tracks) { + std::unique_ptr fittedTrack(GetFittedTrack(&track, fitMetadata, rawEvent, event)); + trackingEvent->AddFittedTrack(std::move(fittedTrack)); + } } diff --git a/AtReconstruction/AtFitter/AtFitter.h b/AtReconstruction/AtFitter/AtFitter.h index 4414388fc..768572e3d 100644 --- a/AtReconstruction/AtFitter/AtFitter.h +++ b/AtReconstruction/AtFitter/AtFitter.h @@ -20,9 +20,12 @@ class TClass; class TMemberInspector; class AtRawEvent; class AtEvent; +class AtFitMetadata; class AtPatternEvent; +class AtRawEvent; +class AtTrackingEvent; -namespace AtFITTER { +namespace EventFit { class AtFitter : public TObject { public: @@ -31,44 +34,26 @@ class AtFitter : public TObject { using TrackMetadatasSet = std::set>; -protected: - // Pointer to the AtPatternEvent to be fitted. - AtPatternEvent *fPatternEvent{nullptr}; - - // Pointer to the AtFitResult object in which store the fit metadata. - AtFitMetadata *fFitMetadata{nullptr}; - - // Pointers to AtRawEvent and AtEvent. In case some specific fitter needs to access information - // in any of those branches. - AtRawEvent *fRawEvent{nullptr}; - AtEvent *fEvent{nullptr}; - - // Compare function that will be used to sort the fit results for a given track. - virtual bool - CompareTrackFitsFunction(const TrackMetadataPtr &trackMetadataA, const TrackMetadataPtr &trackMetadataB) = 0; - public: AtFitter() = default; ~AtFitter() = default; - virtual std::vector> ProcessEvent() = 0; + virtual void FitEvent(AtTrackingEvent *trackingEvent, AtPatternEvent *patternEvent, + AtFitMetadata *fitMetadata = nullptr, AtRawEvent *rawEvent = nullptr, + AtEvent *event = nullptr); virtual void Init() = 0; - // Mandatory to set. - void SetPatternEvent(AtPatternEvent *patternEvent) { fPatternEvent = patternEvent; } - void SetFitMetadata(AtFitMetadata *fitMetadata) { fFitMetadata = fitMetadata; } - - // Optional to set. - void SetRawEvent(AtRawEvent *rawEvent) { fRawEvent = rawEvent; } - void SetEvent(AtEvent *event) { fEvent = event; } +protected: + virtual AtFittedTrack *GetFittedTrack(AtTrack *track, AtFitMetadata *fitMetadata = nullptr, + AtRawEvent *rawEvent = nullptr, AtEvent *event = nullptr) = 0; - // Reset pointers to AtPatternEvent, AtRawEvent and AtEvent. - void Reset(); + // Compare function that will be used to sort the fit results for a given track. + virtual bool + CompareTrackFitsFunction(const TrackMetadataPtr &trackMetadataA, const TrackMetadataPtr &trackMetadataB) = 0; -protected: ClassDef(AtFitter, 2); }; -} // namespace AtFITTER +} // namespace EventFit #endif diff --git a/AtReconstruction/AtFitter/AtGenfit.h b/AtReconstruction/AtFitter/AtGenfit.h index e42aaa271..60c97ba59 100644 --- a/AtReconstruction/AtFitter/AtGenfit.h +++ b/AtReconstruction/AtFitter/AtGenfit.h @@ -1,7 +1,7 @@ #ifndef ATGENFIT_H #define ATGENFIT_H -#include "AtFitter.h" +#include "AtFitterOld.h" #include "AtFormat.h" #include "AtKinematics.h" #include "AtParsers.h" @@ -55,7 +55,9 @@ class MeasurementFactory; namespace AtFITTER { -class AtGenfit : public AtFitter { +class [[deprecated( + "This still derives from the old AtFitter. Please consider updating it to the new AtFitter at some point")]] AtGenfit + : public AtFitterOld { private: std::shared_ptr fKalmanFitter; TClonesArray *fGenfitTrackArray; diff --git a/AtReconstruction/AtFitterTask.cxx b/AtReconstruction/AtFitterTask.cxx index aaf52d044..dd6e4b940 100644 --- a/AtReconstruction/AtFitterTask.cxx +++ b/AtReconstruction/AtFitterTask.cxx @@ -26,7 +26,7 @@ class AtFittedTrack; ClassImp(AtFitterTask); -AtFitterTask::AtFitterTask(std::unique_ptr fitter) +AtFitterTask::AtFitterTask(std::unique_ptr fitter) : fInputBranchName("AtPatternEvent"), fOutputBranchName("AtTrackingEvent"), fIsPersistence(kFALSE), fTrackingEventArray(TClonesArray("AtTrackingEvent", 1)), fFitter(std::move(fitter)), fRawEventBranchName(""), fEventBranchName(""), fFitMetadataBranchName("") @@ -113,42 +113,35 @@ void AtFitterTask::SetParContainers() void AtFitterTask::Exec(Option_t *option) { - fFitter->Reset(); - if (fPatternEventArray->GetEntriesFast() == 0) return; - // If there is AtRawEvent available, pass it to the fitter. - if (fRawEventArray) { - AtRawEvent *rawEvent = dynamic_cast(fRawEventArray->At(0)); - fFitter->SetRawEvent(rawEvent); - } + // If there is AtRawEvent available, get it so it can be passed to the fitter. + AtRawEvent *rawEvent = nullptr; + if (fRawEventArray) + rawEvent = dynamic_cast(fRawEventArray->At(0)); - // If there is AtEvent available, pass it to the fitter. - if (fEventArray) { - AtEvent *event = dynamic_cast(fEventArray->At(0)); - fFitter->SetEvent(event); - } + // If there is AtEvent available, get it so it can be passed to the fitter. + AtEvent *event = nullptr; + if (fEventArray) + event = dynamic_cast(fEventArray->At(0)); fTrackingEventArray.Delete(); auto trackingEvent = dynamic_cast(fTrackingEventArray.ConstructedAt(0)); auto fitMetadata = dynamic_cast(fFitMetadataArray.ConstructedAt(0)); - LOG(info) << " AtFitterTask::Exec() : Fitting event " << fEventCnt; + LOG(info) << " Fitting event " << fEventCnt; AtPatternEvent *patternEvent = dynamic_cast(fPatternEventArray->At(0)); std::vector &tracks = patternEvent->GetTrackCand(); - LOG(info) << " AtFitterTask::Exec() : Number of candidate tracks : " << tracks.size(); + LOG(info) << " Number of candidate tracks : " << tracks.size(); - fFitter->SetPatternEvent(patternEvent); - fFitter->SetFitMetadata(fitMetadata); - auto fittedTracks = fFitter->ProcessEvent(); + fFitter->FitEvent(trackingEvent, patternEvent, fitMetadata, rawEvent, event); - LOG(info) << " AtFitterTask::Exec() : Number of fitted tracks : " << fittedTracks.size(); + auto &fittedTracks = trackingEvent->GetFittedTracks(); - for (auto &fittedTrack : fittedTracks) - trackingEvent->AddFittedTrack(std::move(fittedTrack)); + LOG(info) << " Number of fitted tracks : " << fittedTracks.size(); ++fEventCnt; } diff --git a/AtReconstruction/AtFitterTask.h b/AtReconstruction/AtFitterTask.h index c950382c0..0537cd80d 100644 --- a/AtReconstruction/AtFitterTask.h +++ b/AtReconstruction/AtFitterTask.h @@ -37,14 +37,14 @@ class AtTrack; namespace AtTools { class AtTrackTransformer; } // namespace AtTools -namespace AtFITTER { +namespace EventFit { class AtFitter; -} // namespace AtFITTER +} // namespace EventFit class AtFitterTask : public FairTask { public: - AtFitterTask(std::unique_ptr fitter); + AtFitterTask(std::unique_ptr fitter); ~AtFitterTask() = default; void SetInputBranch(TString branchName); @@ -66,7 +66,7 @@ class AtFitterTask : public FairTask { Bool_t fIsPersistence; //!< Persistence check variable - std::unique_ptr fFitter; + std::unique_ptr fFitter; AtDigiPar *fPar{nullptr}; TClonesArray *fPatternEventArray; TClonesArray fTrackingEventArray; diff --git a/AtReconstruction/AtReconstructionLinkDef.h b/AtReconstruction/AtReconstructionLinkDef.h index bd4887400..384b6611f 100755 --- a/AtReconstruction/AtReconstructionLinkDef.h +++ b/AtReconstruction/AtReconstructionLinkDef.h @@ -46,14 +46,15 @@ #pragma link C++ namespace kf; #pragma link C++ namespace kf::util; -#pragma link C++ namespace AtFITTER; +#pragma link C++ namespace EventFit; #pragma link C++ class AtFitterTask + ; -#pragma link C++ class AtFITTER::AtFitter - ; +#pragma link C++ class EventFit::AtFitter - ; /* Classes that depend on Genfit2 */ +#pragma link C++ namespace AtFITTER; #pragma link C++ class genfit::AtSpacepointMeasurement + ; #pragma link C++ class AtFITTER::AtFitterOld + ; -#pragma link C++ class AtFITTER::AtFitterTaskOld + ; +#pragma link C++ class AtFitterTaskOld + ; #pragma link C++ class AtFITTER::AtGenfit + ; #pragma link C++ namespace MCFitter; From ed58fd9f332fd165c64f6e913495998438c136e0 Mon Sep 17 00:00:00 2001 From: RealAurio Date: Wed, 3 Sep 2025 13:04:59 +0200 Subject: [PATCH 7/9] Last change? --- AtReconstruction/AtFitterTask.h | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/AtReconstruction/AtFitterTask.h b/AtReconstruction/AtFitterTask.h index 0537cd80d..41c1afa7f 100644 --- a/AtReconstruction/AtFitterTask.h +++ b/AtReconstruction/AtFitterTask.h @@ -41,8 +41,13 @@ namespace EventFit { class AtFitter; } // namespace EventFit +/** + * Task that takes a certain AtFitter and uses to fit an AtPatternEvent. The AtFitter may need access to the AtRawEvent + * or AtEvent as well, so pointers to them are also read and passed to the AtFitter. An AtFitMetadata object may also be + * written, which would contain the fit metadata information for all fits done to all AtTracks. + * Specific logic of the fitting is contained in AtFitter and derived classes. + */ class AtFitterTask : public FairTask { - public: AtFitterTask(std::unique_ptr fitter); ~AtFitterTask() = default; @@ -53,7 +58,6 @@ class AtFitterTask : public FairTask { void SetRawEventBranch(TString branchName); void SetEventBranch(TString branchName); - void SetFitMetadataBranch(TString branchName); virtual InitStatus Init(); From 841e510a98ccadca64c94c1b069f692b153e524c Mon Sep 17 00:00:00 2001 From: RealAurio Date: Fri, 5 Sep 2025 16:44:39 +0200 Subject: [PATCH 8/9] Very last change? --- AtReconstruction/AtFitter/AtFitter.cxx | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/AtReconstruction/AtFitter/AtFitter.cxx b/AtReconstruction/AtFitter/AtFitter.cxx index 81e02ef4d..03f0112fc 100644 --- a/AtReconstruction/AtFitter/AtFitter.cxx +++ b/AtReconstruction/AtFitter/AtFitter.cxx @@ -8,6 +8,18 @@ ClassImp(EventFit::AtFitter); void EventFit::AtFitter::FitEvent(AtTrackingEvent *trackingEvent, AtPatternEvent *patternEvent, AtFitMetadata *fitMetadata, AtRawEvent *rawEvent, AtEvent *event) { + // Check for nullptr. + if (trackingEvent == nullptr) { + LOG(error) << " Tracking event is nullptr! The fitter can not fit this event. Maybe the tracking event is not " + "being constructed properly in the fitter task."; + return; + } + + if (patternEvent == nullptr) { + LOG(error) << " Pattern event is nullptr! The fitter can not fit this event."; + return; + } + // Extract the candidate AtTracks. std::vector tracks = patternEvent->GetTrackCand(); From 11b7a6f9d4780ee184c71d08062031ad6efb0fb5 Mon Sep 17 00:00:00 2001 From: RealAurio Date: Mon, 22 Sep 2025 12:56:31 +0200 Subject: [PATCH 9/9] Removed stuff from AtFitter that ended up being too case specific to make sense to keep in the base class. Also fixed some serious bugs that I found while actually implementing a fitter class in the braggCurveFitter feature branch. --- AtData/AtFitMetadata.h | 2 +- AtReconstruction/AtFitter/AtFitter.cxx | 6 +++--- AtReconstruction/AtFitter/AtFitter.h | 12 ------------ AtReconstruction/AtFitterTask.cxx | 5 ++--- AtReconstruction/AtFitterTask.h | 2 -- AtReconstruction/AtReconstructionLinkDef.h | 2 +- 6 files changed, 7 insertions(+), 22 deletions(-) diff --git a/AtData/AtFitMetadata.h b/AtData/AtFitMetadata.h index 9109e0e38..da8fa1c3c 100644 --- a/AtData/AtFitMetadata.h +++ b/AtData/AtFitMetadata.h @@ -45,7 +45,7 @@ class AtFitMetadata : public TObject { void SetTrackMetadatasVector(Int_t trackID, TrackMetadatasVector metadatas) { - fMetadatas[trackID] = std::move(metadatas); + fMetadatas.insert({trackID, std::move(metadatas)}); } void SetEventID(ULong_t id) { fEventID = id; } diff --git a/AtReconstruction/AtFitter/AtFitter.cxx b/AtReconstruction/AtFitter/AtFitter.cxx index 03f0112fc..2dd2582d9 100644 --- a/AtReconstruction/AtFitter/AtFitter.cxx +++ b/AtReconstruction/AtFitter/AtFitter.cxx @@ -3,8 +3,6 @@ #include "AtPatternEvent.h" #include "AtTrackingEvent.h" -ClassImp(EventFit::AtFitter); - void EventFit::AtFitter::FitEvent(AtTrackingEvent *trackingEvent, AtPatternEvent *patternEvent, AtFitMetadata *fitMetadata, AtRawEvent *rawEvent, AtEvent *event) { @@ -20,8 +18,10 @@ void EventFit::AtFitter::FitEvent(AtTrackingEvent *trackingEvent, AtPatternEvent return; } - // Extract the candidate AtTracks. + // Extract the candidate AtTracks. If there are not any tracks, return earlier. std::vector tracks = patternEvent->GetTrackCand(); + if (!tracks.size()) + return; // Save the original AtTracks to the AtTrackingEvent. trackingEvent->SetTrackArray(&tracks); diff --git a/AtReconstruction/AtFitter/AtFitter.h b/AtReconstruction/AtFitter/AtFitter.h index 768572e3d..18b8f5190 100644 --- a/AtReconstruction/AtFitter/AtFitter.h +++ b/AtReconstruction/AtFitter/AtFitter.h @@ -28,12 +28,6 @@ class AtTrackingEvent; namespace EventFit { class AtFitter : public TObject { -public: - using TrackMetadataPtr = std::unique_ptr; - using TrackMetadatasVector = std::vector; - using TrackMetadatasSet = - std::set>; - public: AtFitter() = default; ~AtFitter() = default; @@ -46,12 +40,6 @@ class AtFitter : public TObject { protected: virtual AtFittedTrack *GetFittedTrack(AtTrack *track, AtFitMetadata *fitMetadata = nullptr, AtRawEvent *rawEvent = nullptr, AtEvent *event = nullptr) = 0; - - // Compare function that will be used to sort the fit results for a given track. - virtual bool - CompareTrackFitsFunction(const TrackMetadataPtr &trackMetadataA, const TrackMetadataPtr &trackMetadataB) = 0; - - ClassDef(AtFitter, 2); }; } // namespace EventFit diff --git a/AtReconstruction/AtFitterTask.cxx b/AtReconstruction/AtFitterTask.cxx index dd6e4b940..c4ce51af5 100644 --- a/AtReconstruction/AtFitterTask.cxx +++ b/AtReconstruction/AtFitterTask.cxx @@ -24,12 +24,10 @@ class AtTrack; class AtFittedTrack; -ClassImp(AtFitterTask); - AtFitterTask::AtFitterTask(std::unique_ptr fitter) : fInputBranchName("AtPatternEvent"), fOutputBranchName("AtTrackingEvent"), fIsPersistence(kFALSE), fTrackingEventArray(TClonesArray("AtTrackingEvent", 1)), fFitter(std::move(fitter)), fRawEventBranchName(""), - fEventBranchName(""), fFitMetadataBranchName("") + fEventBranchName(""), fFitMetadataBranchName(""), fFitMetadataArray(TClonesArray("AtFitMetadata", 1)) { } @@ -127,6 +125,7 @@ void AtFitterTask::Exec(Option_t *option) event = dynamic_cast(fEventArray->At(0)); fTrackingEventArray.Delete(); + fFitMetadataArray.Delete(); auto trackingEvent = dynamic_cast(fTrackingEventArray.ConstructedAt(0)); auto fitMetadata = dynamic_cast(fFitMetadataArray.ConstructedAt(0)); diff --git a/AtReconstruction/AtFitterTask.h b/AtReconstruction/AtFitterTask.h index 41c1afa7f..771ad9a8d 100644 --- a/AtReconstruction/AtFitterTask.h +++ b/AtReconstruction/AtFitterTask.h @@ -87,8 +87,6 @@ class AtFitterTask : public FairTask { TString fFitMetadataBranchName; TClonesArray fFitMetadataArray; bool fSaveFitMetadata{false}; - - ClassDef(AtFitterTask, 2); }; #endif diff --git a/AtReconstruction/AtReconstructionLinkDef.h b/AtReconstruction/AtReconstructionLinkDef.h index 384b6611f..8e543d0ac 100755 --- a/AtReconstruction/AtReconstructionLinkDef.h +++ b/AtReconstruction/AtReconstructionLinkDef.h @@ -48,7 +48,7 @@ #pragma link C++ namespace EventFit; #pragma link C++ class AtFitterTask + ; -#pragma link C++ class EventFit::AtFitter - ; +#pragma link C++ class EventFit::AtFitter - !; /* Classes that depend on Genfit2 */ #pragma link C++ namespace AtFITTER;