src/smear/Bremsstrahlung.cxx

0001 /**
0002  \file
0003  Implementation of class Smear::Bremsstrahlung.
0004
0005  \author    Michael Savastio
0006  \date      2011-08-19
0007  \copyright 2011 Brookhaven National Lab
0008  */
0009
0010 #include "eicsmear/smear/Bremsstrahlung.h"
0011
0012 #include <cassert>
0013 #include <cmath>
0014
0015 #include "eicsmear/erhic/VirtualParticle.h"
0016
0017 namespace Smear {
0018
0019 Bremsstrahlung::Bremsstrahlung(double epsilon,
0020                                double traversed,
0021                                double radLength)
0022 : mParticle(nullptr)
0023 , mKMin(0.)
0024 , mKMax(0.)
0025 , mEpsilon(epsilon)
0026 , mTraversed(traversed)
0027 , mRadLength(radLength)
0028 , mPdf(NULL) {
0029   Accept.AddParticle(11);
0030   Accept.AddParticle(-11);
0031 }
0032
0033 Bremsstrahlung::Bremsstrahlung(const Bremsstrahlung& other)
0034     : Device(other), mParticle(nullptr), mKMin(other.mKMin), mKMax(other.mKMax),
0035       mEpsilon(other.mEpsilon), mTraversed(other.mTraversed),
0036       mRadLength(other.mRadLength), mPdf(NULL) {
0037   // Duplicate cached particle if there is one
0038   // SetParticle() duplicates the particle and sets up the PDF
0039   //if (other.mParticle.get()) {
0040   //  SetParticle(*mParticle);
0041   //}  // if
0042 }
0043
0044 double Bremsstrahlung::dSigmadK(double *x, double*) {
0045   double k = x[0];
0046   double ret = 4. / 3.;
0047   ret += -4. * k / (3. * mParticle->E);
0048   ret += pow(k / mParticle->E, 2.);
0049   ret /= k;
0050   return ret;
0051 }
0052
0053 bool Bremsstrahlung::SetupPDF() {
0054   double lower = mEpsilon;
0055   double upper = mParticle->E - mEpsilon;
0056   if (upper < lower || std::isnan(upper) || std::isnan(lower)) {
0057     return false;
0058   }  // if
0059   mKMin = lower;
0060   mKMax = upper;
0061   mPdf->SetRange(mKMin, mKMax);
0062   return true;
0063 }
0064
0065 int Bremsstrahlung::NGamma() {
0066   double ret = 4. * log(mKMax / mKMin) / 3.;
0067   ret += -4. * (mKMax - mKMin) / (3. * mParticle->E);
0068   ret += 0.5* pow((mKMax - mKMin) / mParticle->E, 2.);
0069   ret *= mTraversed / mRadLength;
0070   int n = static_cast<int>(ret);
0071   if (fabs(ret - n) < fabs(ret - n - 1)) {
0072     return n;
0073   } else {
0074     return n + 1;
0075   }  // if
0076 }
0077
0078 void Bremsstrahlung::SetParticle(const erhic::VirtualParticle& prt) {
0079   mParticle.reset(static_cast<erhic::ParticleMC*>(prt.Clone()));
0080   if (!mPdf) {
0081     mPdf = new TF1("Smear_Bremsstrahlung_PDF",
0082                    this,
0083                    &Smear::Bremsstrahlung::dSigmadK,
0084                    mKMin, mKMax, 0);
0085   }  // if
0086   SetupPDF();
0087 }
0088
0089 void Bremsstrahlung::FixParticleKinematics(ParticleMCS& prt) {
0090   prt.SetP(sqrt(prt.GetE() * prt.GetE() - prt.GetM() * prt.GetM()) );
0091   if (prt.GetP() < 0. || std::isnan(prt.GetP())) prt.SetP(0.);
0092   prt.SetPt( prt.GetP() * sin(prt.GetTheta() ));
0093   prt.SetPz( prt.GetP() * cos(prt.GetTheta() ));
0094 }
0095
0096 Bremsstrahlung* Bremsstrahlung::Clone(Option_t* /* not used */) const {
0097   return new Bremsstrahlung(*this);
0098 }
0099
0100 void Bremsstrahlung::Smear(const erhic::VirtualParticle& prt,
0101                            ParticleMCS& prtOut) {
0102   SetParticle(prt);
0103   const int nGamma = NGamma();
0104   for (int i = 0; i < nGamma; i++) {
0105     if (!SetupPDF()) break;
0106     double loss = mPdf->GetRandom();
0107     mParticle->E -= loss;
0108   }  // for
0109   prtOut.SetE ( mParticle->GetE() );
0110   FixParticleKinematics(prtOut);
0111   prtOut.HandleBogusValues(kE);
0112   mParticle.reset(NULL);
0113 }
0114
0115 }  // namespace Smear