Bremsstrahlung.cxx

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Bremsstrahlung.cxx

#include "eicsmear/smear/Bremsstrahlung.h"

#include <cassert>
#include <cmath>

#include "eicsmear/erhic/VirtualParticle.h"

namespace Smear {

Bremsstrahlung::Bremsstrahlung(double epsilon,
                               double traversed,
                               double radLength)
: mParticle(nullptr)
, mKMin(0.)
, mKMax(0.)
, mEpsilon(epsilon)
, mTraversed(traversed)
, mRadLength(radLength)
, mPdf(NULL) {
  Accept.AddParticle(11);
  Accept.AddParticle(-11);
}

Bremsstrahlung::Bremsstrahlung(const Bremsstrahlung& other)
    : Device(other), mParticle(nullptr), mKMin(other.mKMin), mKMax(other.mKMax),
      mEpsilon(other.mEpsilon), mTraversed(other.mTraversed),
      mRadLength(other.mRadLength), mPdf(NULL) {
  // Duplicate cached particle if there is one
  // SetParticle() duplicates the particle and sets up the PDF
  //if (other.mParticle.get()) {
  //  SetParticle(*mParticle);
  //}  // if
}

double Bremsstrahlung::dSigmadK(double *x, double*) {
  double k = x[0];
  double ret = 4. / 3.;
  ret += -4. * k / (3. * mParticle->E);
  ret += pow(k / mParticle->E, 2.);
  ret /= k;
  return ret;
}

bool Bremsstrahlung::SetupPDF() {
  double lower = mEpsilon;
  double upper = mParticle->E - mEpsilon;
  if (upper < lower || std::isnan(upper) || std::isnan(lower)) {
    return false;
  }  // if
  mKMin = lower;
  mKMax = upper;
  mPdf->SetRange(mKMin, mKMax);
  return true;
}

int Bremsstrahlung::NGamma() {
  double ret = 4. * log(mKMax / mKMin) / 3.;
  ret += -4. * (mKMax - mKMin) / (3. * mParticle->E);
  ret += 0.5* pow((mKMax - mKMin) / mParticle->E, 2.);
  ret *= mTraversed / mRadLength;
  int n = static_cast<int>(ret);
  if (fabs(ret - n) < fabs(ret - n - 1)) {
    return n;
  } else {
    return n + 1;
  }  // if
}

void Bremsstrahlung::SetParticle(const erhic::VirtualParticle& prt) {
  mParticle.reset(static_cast<erhic::ParticleMC*>(prt.Clone()));
  if (!mPdf) {
    mPdf = new TF1("Smear_Bremsstrahlung_PDF",
                   this,
                   &Smear::Bremsstrahlung::dSigmadK,
                   mKMin, mKMax, 0);
  }  // if
  SetupPDF();
}

void Bremsstrahlung::FixParticleKinematics(ParticleMCS& prt) {
  prt.SetP(sqrt(prt.GetE() * prt.GetE() - prt.GetM() * prt.GetM()) );
  if (prt.GetP() < 0. || std::isnan(prt.GetP())) prt.SetP(0.);
  prt.SetPt( prt.GetP() * sin(prt.GetTheta() ));
  prt.SetPz( prt.GetP() * cos(prt.GetTheta() ));
}

Bremsstrahlung* Bremsstrahlung::Clone(Option_t* /* not used */) const {
  return new Bremsstrahlung(*this);
}

void Bremsstrahlung::Smear(const erhic::VirtualParticle& prt,
                           ParticleMCS& prtOut) {
  SetParticle(prt);
  const int nGamma = NGamma();
  for (int i = 0; i < nGamma; i++) {
    if (!SetupPDF()) break;
    double loss = mPdf->GetRandom();
    mParticle->E -= loss;
  }  // for
  prtOut.SetE ( mParticle->GetE() );
  FixParticleKinematics(prtOut);
  prtOut.HandleBogusValues(kE);
  mParticle.reset(NULL);
}

}  // namespace Smear