/
Bremsstrahlung.cxx
115 lines (100 loc) · 3.04 KB
/
Bremsstrahlung.cxx
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
/**
\file
Implementation of class Smear::Bremsstrahlung.
\author Michael Savastio
\date 2011-08-19
\copyright 2011 Brookhaven National Lab
*/
#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