Skip to content

Commit

Permalink
Merge pull request #109 from rest-for-physics/missing-observables
Browse files Browse the repository at this point in the history 
TRestRawSignal - Use double for storage instead of float
  • Loading branch information
@lobis
lobis authored Mar 9, 2024
2 parents 3a0131c + c15ca0d commit 75c232b
Show file tree
Hide file tree
Showing 4 changed files with 45 additions and 35 deletions.
16 changes: 8 additions & 8 deletions inc/TRestDetectorSignal.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,8 +37,8 @@ class TRestDetectorSignal {
protected:
Int_t fSignalID = -1;

std::vector<Float_t> fSignalTime; // Vector with the time of the signal
std::vector<Float_t> fSignalCharge; // Vector with the charge of the signal
std::vector<Double_t> fSignalTime; // Vector with the time of the signal
std::vector<Double_t> fSignalCharge; // Vector with the charge of the signal

// TODO: remove this and use readout
std::string fName; // Name of the signal
Expand All @@ -52,7 +52,7 @@ class TRestDetectorSignal {
void IncreaseAmplitude(const TVector2& p);
void SetPoint(const TVector2& p);

// TODO other objects should probably skip using GetMaxIndex direclty
// TODO other objects should probably skip using GetMaxIndex directly
Int_t GetMaxIndex(Int_t from = 0, Int_t to = 0);

TVector2 GetMaxGauss();
Expand Down Expand Up @@ -97,7 +97,7 @@ class TRestDetectorSignal {

void Normalize(Double_t scale = 1.);

std::vector<Int_t> GetPointsOverThreshold() { return fPointsOverThreshold; }
std::vector<Int_t> GetPointsOverThreshold() const { return fPointsOverThreshold; }

Double_t GetAverage(Int_t start = 0, Int_t end = 0);
Int_t GetMaxPeakWidth();
Expand All @@ -114,14 +114,14 @@ class TRestDetectorSignal {
Double_t GetMinTime() const;
Double_t GetMaxTime() const;

Double_t GetData(Int_t index) const { return (double)fSignalCharge[index]; }
Double_t GetTime(Int_t index) const { return (double)fSignalTime[index]; }
Double_t GetData(Int_t index) const { return fSignalCharge[index]; }
Double_t GetTime(Int_t index) const { return fSignalTime[index]; }

// Setters
void SetSignalID(Int_t sID) { fSignalID = sID; }
void SetID(Int_t sID) { fSignalID = sID; }

void NewPoint(Float_t time, Float_t data);
void NewPoint(Double_t time, Double_t data);
void IncreaseAmplitude(Double_t t, Double_t d);

void SetPoint(Double_t t, Double_t d);
Expand Down Expand Up @@ -167,6 +167,6 @@ class TRestDetectorSignal {
// Destructor
~TRestDetectorSignal();

ClassDef(TRestDetectorSignal, 3);
ClassDef(TRestDetectorSignal, 4);
};
#endif
13 changes: 13 additions & 0 deletions src/TRestDetectorHitsToSignalProcess.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -267,6 +267,13 @@ TRestEvent* TRestDetectorHitsToSignalProcess::ProcessEvent(TRestEvent* inputEven
velocity = REST_Physics::lightSpeed;
}

if (velocity <= 0) {
RESTError
<< "TRestDetectorHitsToSignalProcess: Negative velocity. This should not happen."
<< RESTendl;
exit(1);
}

Double_t time = t + distance / velocity;

if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug && hit < 20) {
Expand All @@ -290,6 +297,12 @@ TRestEvent* TRestDetectorHitsToSignalProcess::ProcessEvent(TRestEvent* inputEven

time = floor(time / fSampling) * fSampling;

if (time < 0) {
RESTError << "TRestDetectorHitsToSignalProcess: Negative time. This should not happen. "
"EventID: "
<< fHitsEvent->GetID() << RESTendl;
exit(1);
}
fSignalEvent->AddChargeToSignal(daqId, time, energy);

auto signal = fSignalEvent->GetSignalById(daqId);
Expand Down
48 changes: 22 additions & 26 deletions src/TRestDetectorSignal.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,11 +48,9 @@ TRestDetectorSignal::TRestDetectorSignal() {
fPointsOverThreshold.clear();
}

TRestDetectorSignal::~TRestDetectorSignal() {
// TRestDetectorSignal destructor
}
TRestDetectorSignal::~TRestDetectorSignal() = default;

void TRestDetectorSignal::NewPoint(Float_t time, Float_t data) {
void TRestDetectorSignal::NewPoint(Double_t time, Double_t data) {
fSignalTime.push_back(time);
fSignalCharge.push_back(data);
}
Expand All @@ -61,10 +59,7 @@ void TRestDetectorSignal::NewPoint(Float_t time, Float_t data) {
/// \brief If the point already exists inside the detector signal event,
/// the amplitude value will be added to the corresponding time.
///
void TRestDetectorSignal::IncreaseAmplitude(Double_t t, Double_t d) {
TVector2 p(t, d);
IncreaseAmplitude(p);
}
void TRestDetectorSignal::IncreaseAmplitude(Double_t t, Double_t d) { IncreaseAmplitude({t, d}); }

///////////////////////////////////////////////
/// \brief If the point already exists inside the detector signal event,
Expand All @@ -73,9 +68,9 @@ void TRestDetectorSignal::IncreaseAmplitude(Double_t t, Double_t d) {
/// The input vector should contain a physical time and an amplitude.
///
void TRestDetectorSignal::IncreaseAmplitude(const TVector2& p) {
Int_t index = GetTimeIndex(p.X());
Float_t x = p.X();
Float_t y = p.Y();
Double_t x = p.X();
Double_t y = p.Y();
Int_t index = GetTimeIndex(x);

if (index >= 0) {
fSignalTime[index] = x;
Expand All @@ -89,14 +84,14 @@ void TRestDetectorSignal::IncreaseAmplitude(const TVector2& p) {
///////////////////////////////////////////////
/// \brief If the point already exists inside the detector signal event,
/// it will be overwritten. If it does not exists, a new point will be
/// added to the poins vector.
/// added to the points vector.
///
/// The input vector should contain a physical time and an amplitude.
///
void TRestDetectorSignal::SetPoint(const TVector2& p) {
Int_t index = GetTimeIndex(p.X());
Float_t x = p.X();
Float_t y = p.Y();
Double_t x = p.X();
Double_t y = p.Y();

if (index >= 0) {
fSignalTime[index] = x;
Expand Down Expand Up @@ -471,7 +466,7 @@ Double_t TRestDetectorSignal::GetMaxPeakTime(Int_t from, Int_t to) { return GetT
Double_t TRestDetectorSignal::GetMinPeakValue() { return GetData(GetMinIndex()); }

Int_t TRestDetectorSignal::GetMinIndex() const {
Double_t min = numeric_limits<double>::max();
Double_t min = numeric_limits<Double_t>::max();
Int_t index = 0;

for (int i = 0; i < GetNumberOfPoints(); i++) {
Expand All @@ -485,35 +480,35 @@ Int_t TRestDetectorSignal::GetMinIndex() const {
}

Double_t TRestDetectorSignal::GetMinTime() const {
Double_t minTime = numeric_limits<float>::max();
if (GetNumberOfPoints() == 0) {
return 0;
}
Double_t minTime = numeric_limits<Double_t>::max();
for (int i = 0; i < GetNumberOfPoints(); i++) {
const auto time = GetTime(i);
const Double_t time = GetTime(i);
if (time < minTime) {
minTime = time;
}
}
if (minTime == numeric_limits<float>::max()) {
minTime = 0;
}
return minTime;
}

Double_t TRestDetectorSignal::GetMaxTime() const {
Double_t maxTime = numeric_limits<float>::min();
if (GetNumberOfPoints() == 0) {
return 0;
}
Double_t maxTime = numeric_limits<Double_t>::min();
for (int i = 0; i < GetNumberOfPoints(); i++) {
const auto time = GetTime(i);
if (time > maxTime) {
maxTime = time;
}
}
if (maxTime == numeric_limits<float>::min()) {
maxTime = 0;
}
return maxTime;
}

Int_t TRestDetectorSignal::GetTimeIndex(Double_t t) {
Float_t time = t;
Double_t time = t;

for (int n = 0; n < GetNumberOfPoints(); n++) {
if (time == fSignalTime[n]) {
Expand Down Expand Up @@ -635,9 +630,10 @@ void TRestDetectorSignal::MultiplySignalBy(Double_t factor) {

void TRestDetectorSignal::ExponentialConvolution(Double_t fromTime, Double_t decayTime, Double_t offset) {
for (int i = 0; i < GetNumberOfPoints(); i++) {
if (fSignalTime[i] > fromTime)
if (fSignalTime[i] > fromTime) {
fSignalCharge[i] =
(fSignalCharge[i] - offset) * exp(-(fSignalTime[i] - fromTime) / decayTime) + offset;
}
}
}

Expand Down
3 changes: 2 additions & 1 deletion src/TRestDetectorSignalEvent.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -138,10 +138,11 @@ Double_t TRestDetectorSignalEvent::GetMinValue() {

Double_t TRestDetectorSignalEvent::GetMinTime() {
Double_t minTime = numeric_limits<Double_t>::max();
for (int s = 0; s < GetNumberOfSignals(); s++)
for (int s = 0; s < GetNumberOfSignals(); s++) {
if (minTime > fSignal[s].GetMinTime()) {
minTime = fSignal[s].GetMinTime();
}
}
return minTime;
}

Expand Down

0 comments on commit 75c232b

Please sign in to comment.