Fit speed develop

.gitignore

@@ -60,4 +60,5 @@ macros/data
 compiled/E20009Analysis/.NabinWorkDir/
 
 .claude/plans/*
+.cache
 

AtAnalysis/E12014/AtE12014.cxx

@@ -13,7 +13,6 @@
 
 #include <FairLogger.h>
 
-#include <TF1.h>
 #include <TH1.h>
 #include <TString.h>
 
@@ -128,8 +127,8 @@ int E12014::FillHitSums(std::vector<double> &exp, std::vector<double> &sim, cons
       if (fMap->GetPadSize(expHit->GetPadNum()) != 0)
          continue;
 
-      auto funcExp = AtTools::GetHitFunctionTB(*expHit, fPar);
-      if (funcExp == nullptr)
+      auto expGauss = AtTools::GetHitGaussianTB(*expHit, fPar);
+      if (!expGauss)
          continue;
 
       // We have a hit we want to save for both exp and fisison.
@@ -143,8 +142,8 @@ int E12014::FillHitSums(std::vector<double> &exp, std::vector<double> &sim, cons
       }
       if (simHit == nullptr)
          continue;
-      auto funcSim = AtTools::GetHitFunctionTB(*simHit, fPar);
-      if (funcSim == nullptr)
+      auto simGauss = AtTools::GetHitGaussianTB(*simHit, fPar);
+      if (!simGauss)
          continue;
 
       AtPad *pad = nullptr;
@@ -156,10 +155,10 @@ int E12014::FillHitSums(std::vector<double> &exp, std::vector<double> &sim, cons
       //  We now have the sim and exp hits. Fill the arrays
       for (int tb = fTBMin; tb < 512; ++tb) {
 
-         auto val = funcExp->Eval(tb);
+         auto val = (*expGauss)(tb);
          if (val > threshold) {
             exp[tb] += val;
-            sim[tb] += funcSim->Eval(tb);
+            sim[tb] += (*simGauss)(tb);
             if (pad && expADC)
                (*expADC)[tb] += pad->GetADC(tb);
          }
@@ -179,12 +178,12 @@ void E12014::FillSimSum(std::vector<double> &sim, const std::vector<AtHit *> &si
 
    for (auto &expHit : simHits) {
 
-      auto funcExp = AtTools::GetHitFunctionTB(*expHit);
-      if (funcExp == nullptr)
+      auto expGauss = AtTools::GetHitGaussianTB(*expHit);
+      if (!expGauss)
          continue;
 
       for (int tb = fTBMin; tb < 512; ++tb) {
-         sim[tb] += funcExp->Eval(tb);
+         sim[tb] += (*expGauss)(tb);
       }
    }
 }
@@ -269,13 +268,13 @@ E12014::FillHitSum(std::vector<double> &vec, const std::vector<AtHit *> &hits, i
       if (hit->GetCharge() > satThresh)
          continue;
 
-      auto func = AtTools::GetHitFunctionTB(*hit);
-      if (func == nullptr)
+      auto gauss = AtTools::GetHitGaussianTB(*hit);
+      if (!gauss)
          continue;
 
       // Add the charge to the array
       for (int tb = 0; tb < vec.size(); ++tb) {
-         auto val = func->Eval(tb);
+         auto val = (*gauss)(tb);
          if (val > threshold)
             vec[tb] += val;
       }
@@ -309,15 +308,15 @@ void E12014::FillSimHitSum(std::vector<double> &vec, const std::vector<AtHit *>
       if (hit->GetCharge() > satThresh)
          continue;
 
-      auto func = AtTools::GetHitFunctionTB(*hit);
-      if (func == nullptr)
+      auto gauss = AtTools::GetHitGaussianTB(*hit);
+      if (!gauss)
          continue;
 
       // Add the charge to the array
       LOG(debug) << "Adding pad " << hit->GetPadNum();
 
       for (int tb = 0; tb < vec.size(); ++tb) {
-         auto val = func->Eval(tb) * amp;
+         auto val = amp * (*gauss)(tb);
          if (val > threshold)
             vec[tb] += val;
       }

AtData/AtContainerManip.h

@@ -3,7 +3,8 @@
 
 #include <FairLogger.h>
 
-#include <TH1.h> // for TH1, TH1D (ptr only)
+#include <TDirectory.h> // for TDirectory::TContext
+#include <TH1.h>        // for TH1, TH1D (ptr only)
 
 #include <algorithm> // IWYU pragma: keep
 #include <iterator>  // for make_move_iterator
@@ -44,12 +45,12 @@ std::unique_ptr<Hist> CreateHistFromData(const std::string &name, const T &data,
    if (xMin == 0 && xMax == 0)
       xMax = data.size();
 
-   TH1::AddDirectory(false); // Make sure pointer has ownership
+   // Make sure the pointer has ownership. TH1::AddDirectory is a process-wide flag, so don't toggle it from threads.
+   // Build with no current directory instead (gDirectory is thread-local), so the histogram is never appended to one.
+   TDirectory::TContext noDirectory(nullptr);
    auto ret = std::make_unique<Hist>(name.data(), name.data(), data.size(), 0, data.size());
    for (int i = 0; i < ret->GetNbinsX(); ++i)
       ret->SetBinContent(i + 1, data.at(i));
-   // ret->SetDirectory(nullptr);
-   TH1::AddDirectory(true);
    return ret;
 }
 

AtData/AtDataLinkDef.h

@@ -31,6 +31,9 @@
 #pragma link C++ class AtFittedTrack + ;
 #pragma link C++ class AtFissionEvent + ;
 #pragma link C++ class AtGenericTrace + ;
+#pragma link C++ class AtICEvent + ;
+#pragma link C++ class AtENCourseEvent + ;
+#pragma link C++ class AtPPACPair + ;
 
 #pragma link C++ class AtPatterns::AtPattern + ;
 #pragma link C++ class AtPatterns::AtPatternLine + ;

AtData/AtENCourseEvent.cxx

@@ -0,0 +1,38 @@
+#include "AtENCourseEvent.h"
+
+#include <Rtypes.h>
+
+ClassImp(AtENCourseEvent);
+
+AtENCourseEvent::AtENCourseEvent() : AtBaseEvent("AtENCourseEvent")
+{
+   fF2PPACs = nullptr;
+   fF3PPACs = nullptr;
+   for (int i = 0; i < 4; i++)
+      fRFToF[i] = -999;
+   fTDCRef = 0;
+}
+
+AtENCourseEvent::AtENCourseEvent(const AtENCourseEvent &copy) : AtBaseEvent(copy), fTDCRef(copy.fTDCRef)
+{
+   fF2PPACs = copy.fF2PPACs->Clone();
+   fF3PPACs = copy.fF3PPACs->Clone();
+   for (int i = 0; i < 4; i++)
+      fRFToF[i] = copy.fRFToF[i];
+}
+
+AtENCourseEvent &AtENCourseEvent::operator=(AtENCourseEvent object)
+{
+   swap(*this, object);
+   return *this;
+}
+
+void AtENCourseEvent::Clear(Option_t *opt)
+{
+   AtBaseEvent::Clear(opt);
+   fF2PPACs.reset();
+   fF3PPACs.reset();
+   for (int i = 0; i < 4; i++)
+      fRFToF[i] = -999;
+   fTDCRef = 0;
+}

AtData/AtENCourseEvent.h

@@ -0,0 +1,53 @@
+#ifndef AtENCOURSEEVENT_H
+#define AtENCOURSEEVENT_H
+
+#include "AtBaseEvent.h"
+#include "AtPPACPair.h"
+
+#include <FairLogger.h>
+
+#include <Rtypes.h>
+
+class AtENCourseEvent : public AtBaseEvent {
+private:
+   std::unique_ptr<AtPPACPair> fF2PPACs;
+   std::unique_ptr<AtPPACPair> fF3PPACs;
+
+   Int_t fRFToF[4];
+   Int_t fTDCRef;
+
+public:
+   AtENCourseEvent();
+   AtENCourseEvent(const AtENCourseEvent &copy);
+   AtENCourseEvent(const AtBaseEvent &copy) : AtBaseEvent(copy) { SetName("AtENCourseEvent"); }
+   AtENCourseEvent &operator=(const AtENCourseEvent object);
+   virtual ~AtENCourseEvent() = default;
+
+   friend void swap(AtENCourseEvent &first, AtENCourseEvent &second)
+   {
+      using std::swap;
+      swap(dynamic_cast<AtBaseEvent &>(first), dynamic_cast<AtBaseEvent &>(second));
+
+      swap(first.fF2PPACs, second.fF2PPACs);
+      swap(first.fF3PPACs, second.fF3PPACs);
+      swap(first.fRFToF, second.fRFToF);
+      swap(first.fTDCRef, second.fTDCRef);
+   }
+
+   void Clear(Option_t *opt = nullptr) override;
+
+   void SetF2PPACs(std::unique_ptr<AtPPACPair> F2PPACs) { fF2PPACs = std::move(F2PPACs); }
+   void SetF3PPACs(std::unique_ptr<AtPPACPair> F3PPACs) { fF3PPACs = std::move(F3PPACs); }
+   void SetRFToF(Int_t RFToF, Int_t idx = 0) { fRFToF[idx] = RFToF; }
+   void SetTDCRef(Int_t TDCRef) { fTDCRef = TDCRef; }
+
+   const AtPPACPair &GetF2PPACs() const { return *fF2PPACs; }
+   const AtPPACPair &GetF3PPACs() const { return *fF3PPACs; }
+   const Int_t GetRFToF(Int_t idx = 0) const { return fRFToF[idx]; }
+   const Int_t GetTDCRef() const { return fTDCRef; }
+   const Int_t GetTDCRefRFToF(Int_t idx = 0) const { return fRFToF[idx] - fTDCRef; }
+
+   ClassDefOverride(AtENCourseEvent, 1);
+};
+
+#endif

AtData/AtHit.h

@@ -37,6 +37,7 @@ class AtHit : public TObject {
    XYZVector fPositionVariance{}; //< Position variance (unused by AtHitCluster)
    Int_t fHitID;                  //< Unique ID of hit
    Int_t fPadNum;                 //< Pad that generated hit
+   Double_t fCalibratedELoss{-1}; //< Calibrated ELoss from calibration function in AtMap
 
    Double_t fTraceIntegral{-1};     //< Integrated pulse charge of full trace
    Int_t fHitMult{1};               //< Hit multiplicity in the pad where the hit was found
@@ -66,6 +67,7 @@ class AtHit : public TObject {
    virtual void SetPositionVariance(const XYZVector &vec) { fPositionVariance = vec; }
    void SetHitID(Int_t hitID) { fHitID = hitID; }
    void SetPadNum(Int_t padNum) { fPadNum = padNum; }
+   void SetCalibratedELoss(Double_t calibratedELoss) { fCalibratedELoss = calibratedELoss; }
 
    void SetTraceIntegral(Double_t integral) { fTraceIntegral = integral; }
    void SetHitMult(Int_t HitMult) { fHitMult = HitMult; }
@@ -81,6 +83,7 @@ class AtHit : public TObject {
    XYZVector GetPositionSigma() const;
    Double_t GetCharge() const { return fCharge; }
    Int_t GetPadNum() const { return fPadNum; }
+   Double_t GetCalibratedELoss() const { return fCalibratedELoss; }
    Double_t GetTraceIntegral() const { return fTraceIntegral; }
    Int_t GetHitMult() const { return fHitMult; }
    Int_t GetTimeStamp() const { return fTimeStamp; }
@@ -116,7 +119,7 @@ class AtHit : public TObject {
       }
    };
 
-   ClassDef(AtHit, 5);
+   ClassDef(AtHit, 6);
 };
 
 /**

AtData/AtICEvent.cxx

@@ -0,0 +1,17 @@
+#include "AtICEvent.h"
+
+ClassImp(AtICEvent)
+
+   AtICEvent::AtICEvent()
+   : AtBaseEvent("AtICEvent")
+{
+}
+
+void AtICEvent::Clear(Option_t *opt)
+{
+   AtBaseEvent::Clear(opt);
+
+   fRawADC = -1.;
+   fADC = -1.;
+   fIsPedestalSubtracted = kFALSE;
+}

AtData/AtICEvent.h

@@ -0,0 +1,41 @@
+#ifndef AtICEVENT_H
+#define AtICEVENT_H
+
+#include "AtBaseEvent.h"
+
+#include <Rtypes.h>
+
+#include <memory>
+#include <type_traits>
+#include <utility>
+#include <vector>
+class TBuffer;
+class TClass;
+class TMemberInspector;
+
+class AtICEvent : public AtBaseEvent {
+private:
+   Double_t fRawADC = -1.;
+   Double_t fADC = -1.;
+   Bool_t fIsPedestalSubtracted = kFALSE;
+
+public:
+   AtICEvent();
+   AtICEvent(const AtICEvent &) = default;
+   AtICEvent &operator=(const AtICEvent &) = default;
+   virtual ~AtICEvent() = default;
+
+   void Clear(Option_t *opt = nullptr) override;
+
+   void SetRawADC(Double_t value) { fRawADC = value; }
+   void SetADC(Double_t value) { fADC = value; }
+   void SetPedestalSubtracted(Bool_t val = kTRUE) { fIsPedestalSubtracted = val; }
+
+   Double_t GetRawADC() const { return fRawADC; }
+   Double_t GetADC() const { return fADC; }
+   Bool_t IsPedestalSubtracted() const { return fIsPedestalSubtracted; }
+
+   ClassDefOverride(AtICEvent, 1);
+};
+
+#endif

AtData/AtPPACPair.cxx

@@ -0,0 +1,22 @@
+#include "AtPPACPair.h"
+
+ClassImp(AtPPACPair);
+
+void AtPPACPair::CalculateDirection()
+{
+   // Entrance and exit positions relative to the center of the upstream PPAC.
+   XYZPoint initialPoint(fEntrancePosition.X(), fEntrancePosition.Y(), 0);
+   XYZPoint finalPoint(fExitPosition.X(), fExitPosition.Y(), fSeparationDistance);
+
+   // Direction vector.
+   fTrackDirection = finalPoint - initialPoint;
+
+   // Set angles based off the direction vector.
+   fTrackPolarAngle = fTrackDirection.Theta();
+   fTrackAzimutalAngle = fTrackDirection.Phi();
+}
+
+std::unique_ptr<AtPPACPair> AtPPACPair::Clone()
+{
+   return std::make_unique<AtPPACPair>(*this);
+}