QIEInputPulse.cxx

Go to the documentation of this file.

 
#include "TrigScint/QIEInputPulse.h"
 
#include <cmath>
 
namespace trigscint {
 
void QIEInputPulse::addPulse(float toff, float ampl) {
  toff_.push_back(toff);
  ampl_.push_back(ampl);
}
 
float QIEInputPulse::eval(float T) {
  if (ampl_.size() == 0) return 0;
  float val = 0;
  for (int i = 0; i < ampl_.size(); i++) {
    val += evalSingle(T, i);
  }
  return val;
}
 
// bimoid pulse, made out of difference of two sigmoids parametrized by
// rt,ft respectively.
// Parameters:
// rise = rise time
// fall = fall time
Bimoid::Bimoid(float rise, float fall) {
  rt_ = rise;
  ft_ = fall;
}
 
float Bimoid::evalSingle(float T, int id) {
  if (T < toff_[id]) return 0;
  // Normalization constant
  float nc = (ft_ - rt_) * log(2) / ampl_[id];
 
  float y1 = 1 / (1 + exp((toff_[id] - T) / rt_));
  float y2 = 1 / (1 + exp((toff_[id] - T) / ft_));
  return ((y1 - y2) / nc);
}
 
float Bimoid::integrate(float T1, float T2) {
  float val = 0;
  for (int id = 0; id < ampl_.size(); id++) {
    if (ampl_[id] > 0 && T2 > toff_[id]) {
      val += iInt(T2, id) - iInt(T1, id);
    }
  }
  return val;
}
 
float Bimoid::iInt(float T, int id) {
  if (T <= toff_[id]) return 0;
  // Normalization constant
  float nc = (ft_ - rt_) * log(2) / ampl_[id];
 
  float t = T - toff_[id];  // time relative to offset
 
  float ii =  // Integral
      rt_ * log(1 + exp((t - toff_[id]) / rt_)) -
      ft_ * log(1 + exp((t - toff_[id]) / ft_));
 
  return ii / nc;
}
 
float Bimoid::max(int id) {
  float a = 0;
  float b = 50;
  float mx = (a + b) / 2;  // maximum
 
  while (std::abs(derivative(mx, id)) >= 1e-5) {
    if (derivative(a, id) * derivative(mx, id) > 0) {
      a = mx;
    } else
      b = mx;
    mx = (a + b) / 2;
  }
  return (mx);
}
 
float Bimoid::derivative(float T, int id) {
  // Normalization constant
  float nc = (ft_ - rt_) * log(2) / ampl_[id];
 
  float t = T - toff_[id];
  float e1 = exp(-t / rt_);
  float e2 = exp(-t / ft_);
 
  float v1 = e1 / (rt_ * pow(1 + e1, 2));
  float v2 = e2 / (ft_ * pow(1 + e2, 2));
 
  return ((v1 - v2) / nc);  // Actual derivative
}
 
Expo::Expo() {}
 
// A current pulse formed by assuming SiPM as an ideal capacitor which is
// fed with a constant current.
// Parameters:
// k_ = 1/(RC time constant of the capacitor)
// tmax_ = The charging time of the capacitor
Expo::Expo(float k, float tmax) {
  k_ = k;
  tmax_ = tmax;
 
  rt_ = (log(9 + exp(-k_ * tmax_)) - log(1 + 9 * exp(-k_ * tmax_))) / k_;
  ft_ = log(9) / k_;
}
 
// Manually set the rise time and fall time of the pulse
void Expo::setRiseFall(float rr, float ff) {
  rt_ = rr;
  ft_ = ff;
 
  k_ = log(9) / ft_;
  tmax_ = (log(9 - exp(-k_ * rt_)) - log(9 * exp(-k_ * rt_) - 1)) / k_;
}
 
float Expo::evalSingle(float t_, int id) {
  if (id >= ampl_.size()) return 0;
  if (t_ <= toff_[id]) return 0;
  if (ampl_[id] == 0) return 0;
 
  // Normalization constant
  float nc = ampl_[id] / tmax_;
  // time relative to the offset
  float t = t_ - toff_[id];
  if (t < tmax_) {
    return (nc * (1 - exp(-k_ * t)));
  } else {
    return (nc * (1 - exp(-k_ * tmax_)) * exp(k_ * (tmax_ - t)));
  }
  return -1;
}
 
float Expo::max(int id) {
  // Normalization constant
  float nc = ampl_[id] / tmax_;
  return nc * (1 - exp(-k_ * tmax_));
}
 
float Expo::integrate(float T1, float T2) {
  float val = 0;
  for (int id = 0; id < ampl_.size(); id++) {
    if (ampl_[id] > 0 && T2 > toff_[id]) {
      val += iInt(T2, id) - iInt(T1, id);
    }
  }
  return val;
}
 
float Expo::derivative(float T, int id) {
  if (id >= ampl_.size()) return 0;
  if (T <= toff_[id]) return 0;
 
  float t = T - toff_[id];
  // Normalization constant
  float nc = ampl_[id] / tmax_;
 
  if (t <= tmax_) return (nc * k_ * exp(-k_ * t));
  return (-nc * k_ * (1 - exp(-k_ * tmax_)) * exp(k_ * (tmax_ - t)));
}
 
float Expo::iInt(float T, int id) {
  if (T <= toff_[id]) return 0;
  float t = T - toff_[id];
  // Normalization constant
  float nc = ampl_[id] / tmax_;
  if (t < tmax_) return (nc * (k_ * t + exp(-k_ * t) - 1) / k_);
 
  float c1 = (1 - exp(-k_ * tmax_)) / k_;
  float c2 = tmax_ - c1 * exp(k_ * (tmax_ - t));
  return nc * c2;
}
 
}  // namespace trigscint