Skip to content
Snippets Groups Projects
Verified Commit de93a01e authored by Tamas Gal's avatar Tamas Gal :speech_balloon:
Browse files

Cleanup

parent 4595c8f2
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 157 additions and 161 deletions
src/showers.jl
+ 157
− 161
View file @ de93a01e
......@@ -14,48 +14,48 @@ Probability density function for direct light from EM-shower. Returns
- `Δt`: time difference relative to direct Cherenkov light
"""
function directlightfromEMshower(params::LMParameters, pmt::PMTModel, D, cd, θ, ϕ, Δt)
const double ct0 = cd;
const double st0 = sqrt((1.0 + ct0) * (1.0 - ct0));
ct0 = cd;
st0 = sqrt((1.0 + ct0) * (1.0 - ct0));
const double D = std::max(D, getRmin());
const double t = D * getIndexOfRefraction() / C + t_ns; // time [ns]
const double A = getPhotocathodeArea();
D = max(D, getRmin());
t = D * getIndexOfRefraction() / C + t_ns # time [ns]
A = getPhotocathodeArea();
const double px = sin(theta) * cos(phi);
const double pz = cos(theta);
px = sin(theta) * cos(phi);
pz = cos(theta);
const double n0 = getIndexOfRefractionGroup(wmax);
const double n1 = getIndexOfRefractionGroup(wmin);
const double ng = C * t / D; // index of refraction
n0 = getIndexOfRefractionGroup(wmax);
n1 = getIndexOfRefractionGroup(wmin);
ng = C * t / D # index of refraction
if (n0 >= ng) {
return 0.0;
}
end
if (n1 <= ng) {
return 0.0;
}
end
const double w = getWavelength(ng);
const double n = getIndexOfRefractionPhase(w);
w = getWavelength(ng);
n = getIndexOfRefractionPhase(w);
const double l_abs = getAbsorptionLength(w);
const double ls = getScatteringLength(w);
l_abs = getAbsorptionLength(w);
ls = getScatteringLength(w);
const double npe = cherenkov(w, n) * getQE(w);
npe = cherenkov(w, n) * getQE(w);
const double ct = st0 * px + ct0 * pz; // cosine angle of incidence on PMT
ct = st0 * px + ct0 * pz # cosine angle of incidence on PMT
const double U = getAngularAcceptance(ct); // PMT angular acceptance
const double V = exp(-D / l_abs - D / ls); // absorption & scattering
const double W = A / (D * D); // solid angle
U = getAngularAcceptance(ct) # PMT angular acceptance
V = exp(-D / l_abs - D / ls) # absorption & scattering
W = A / (D * D) # solid angle
const double ngp = getDispersionGroup(w);
ngp = getDispersionGroup(w);
const double Ja = D * ngp / C; // dt/dlambda
const double Jb = geant(n, ct0); // d^2N/dcos/dphi
Ja = D * ngp / C # dt/dlambda
Jb = geant(n, ct0) # d^2N/dcos/dphi
return npe * geanc() * U * V * W * Jb / fabs(Ja);
}
end
"""
scatteredlightfromEMshower(params::LMParameters, pmt::PMTModel, D, cd, θ, ϕ, Δt)
......@@ -75,133 +75,133 @@ Probability density function for scattered light from EM-shower. Returns
function scatteredlightfromEMshower(params::LMParameters, pmt::PMTModel, D, cd, θ, ϕ, Δt)
double value = 0;
const double sd = sqrt((1.0 + cd) * (1.0 - cd));
const double D = std::max(D_m, getRmin());
const double L = D;
const double t = D * getIndexOfRefraction() / C + t_ns; // time [ns]
sd = sqrt((1.0 + cd) * (1.0 - cd));
D = max(D_m, getRmin());
L = D;
t = D * getIndexOfRefraction() / C + t_ns # time [ns]
const double A = getPhotocathodeArea();
A = getPhotocathodeArea();
const double px = sin(theta) * cos(phi);
const double py = sin(theta) * sin(phi);
const double pz = cos(theta);
px = sin(theta) * cos(phi);
py = sin(theta) * sin(phi);
pz = cos(theta);
const double qx = cd * px + 0 - sd * pz;
const double qy = 1 * py;
const double qz = sd * px + 0 + cd * pz;
qx = cd * px + 0 - sd * pz;
qy = 1 * py;
qz = sd * px + 0 + cd * pz;
const double n0 = getIndexOfRefractionGroup(wmax);
const double n1 = getIndexOfRefractionGroup(wmin);
n0 = getIndexOfRefractionGroup(wmax);
n1 = getIndexOfRefractionGroup(wmin);
const double ni = C * t / L; // maximal index of refraction
ni = C * t / L # maximal index of refraction
if (n0 >= ni) {
return value;
}
end
const double nj = std::min(ni, n1);
nj = min(ni, n1);
double w = wmax;
w = wmax;
for (const_iterator i = begin(); i != end(); ++i) {
const double ng = 0.5 * (nj + n0) + i->getX() * 0.5 * (nj - n0);
const double dn = i->getY() * 0.5 * (nj - n0);
ng = 0.5 * (nj + n0) + i->getX() * 0.5 * (nj - n0);
dn = i->getY() * 0.5 * (nj - n0);
w = getWavelength(ng, w, 1.0e-5);
const double dw = dn / fabs(getDispersionGroup(w));
dw = dn / fabs(getDispersionGroup(w));
const double n = getIndexOfRefractionPhase(w);
n = getIndexOfRefractionPhase(w);
const double l_abs = getAbsorptionLength(w);
const double ls = getScatteringLength(w);
l_abs = getAbsorptionLength(w);
ls = getScatteringLength(w);
const double npe = cherenkov(w, n) * dw * getQE(w);
npe = cherenkov(w, n) * dw * getQE(w);
if (npe <= 0) {
continue;
}
end
const double Jc = 1.0 / ls; // dN/dx
Jc = 1.0 / ls # dN/dx
const double d = C * t / ng; // photon path
d = C * t / ng # photon path
// const double V = exp(-d/l_abs); //
// V = exp(-d/l_abs) #
// absorption
const double ds = 2.0 / (size() + 1);
ds = 2.0 / (size() + 1);
for (double sb = 0.5 * ds; sb < 1.0 - 0.25 * ds; sb += ds) {
for (int buffer[] = {-1, +1, 0}, *k = buffer; *k != 0; ++k) {
for k in (-1, 1)
const double cb = (*k) * sqrt((1.0 + sb) * (1.0 - sb));
const double dcb = (*k) * ds * sb / cb;
cb = k * sqrt((1.0 + sb) * (1.0 - sb));
dcb = k * ds * sb / cb;
const double v = 0.5 * (d + L) * (d - L) / (d - L * cb);
const double u = d - v;
v = 0.5 * (d + L) * (d - L) / (d - L * cb);
u = d - v;
if (u <= 0) {
continue;
}
end
if (v <= 0) {
continue;
}
end
const double cts = (L * cb - v) / u; // cosine scattering angle
cts = (L * cb - v) / u # cosine scattering angle
const double V =
V =
exp(-d * getInverseAttenuationLength(l_abs, ls, cts));
if (cts < 0.0 &&
v * sqrt((1.0 + cts) * (1.0 - cts)) < MODULE_RADIUS_M) {
continue;
}
end
const double W = std::min(A / (v * v), 2.0 * PI); // solid angle
const double Ja = getScatteringProbability(cts); // d^2P/dcos/dphi
const double Jd = ng * (1.0 - cts) / C; // dt/du
W = min(A / (v * v), 2.0 * PI) # solid angle
Ja = getScatteringProbability(cts) # d^2P/dcos/dphi
Jd = ng * (1.0 - cts) / C # dt/du
const double ca = (L - v * cb) / u;
const double sa = v * sb / u;
ca = (L - v * cb) / u;
sa = v * sb / u;
const double dp = PI / phd.size();
const double dom = dcb * dp * v * v / (u * u);
const double dos = sqrt(dom);
dp = PI / phd.size();
dom = dcb * dp * v * v / (u * u);
dos = sqrt(dom);
for (const_iterator l = phd.begin(); l != phd.end(); ++l) {
const double cp = l->getX();
const double sp = l->getY();
cp = l->getX();
sp = l->getY();
const double ct0 = cd * ca - sd * sa * cp;
ct0 = cd * ca - sd * sa * cp;
const double vx = -sb * cp * qx;
const double vy = -sb * sp * qy;
const double vz = cb * qz;
vx = -sb * cp * qx;
vy = -sb * sp * qy;
vz = cb * qz;
const double U =
U =
getAngularAcceptance(vx + vy + vz) +
getAngularAcceptance(vx - vy + vz); // PMT angular acceptance
getAngularAcceptance(vx - vy + vz) # PMT angular acceptance
// const double Jb = geant(n,ct0); //
// Jb = geant(n,ct0) #
// d^2N/dcos/dphi
// value += npe * geanc() * dom * U * V * W * Ja * Jb * Jc /
// fabs(Jd);
const double Jb = geant(n, ct0 - 0.5 * dos,
ct0 + 0.5 * dos); // dN/dphi
Jb = geant(n, ct0 - 0.5 * dos,
ct0 + 0.5 * dos) # dN/dphi
value += npe * geanc() * dos * U * V * W * Ja * Jb * Jc / fabs(Jd);
}
}
}
}
end
end
end
end
return value;
}
end
"""
directlightfromEMshower(params::LMParameters, pmt::PMTModel, E, D, cd, θ, ϕ, Δt)
......@@ -222,98 +222,98 @@ Probability density function for direct light from EM-shower. Returns
function directlightfromEMshower(params::LMParameters, pmt::PMTModel, E, D, cd, θ, ϕ, Δt)
double value = 0;
const double sd = sqrt((1.0 + cd) * (1.0 - cd));
const double D = std::max(D_m, getRmin());
const double R = D * sd; // minimal distance of approach [m]
const double Z = -D * cd;
const double t = D * getIndexOfRefraction() / C + t_ns; // time [ns]
sd = sqrt((1.0 + cd) * (1.0 - cd));
D = max(D_m, getRmin());
R = D * sd # minimal distance of approach [m]
Z = -D * cd;
t = D * getIndexOfRefraction() / C + t_ns # time [ns]
const double n0 = getIndexOfRefractionGroup(wmax);
const double n1 = getIndexOfRefractionGroup(wmin);
n0 = getIndexOfRefractionGroup(wmax);
n1 = getIndexOfRefractionGroup(wmin);
double zmin = 0.0; // minimal shower length [m]
double zmax = geanz.getLength(E, 1.0); // maximal shower length [m]
double zmin = 0.0 # minimal shower length [m]
double zmax = geanz.getLength(E, 1.0) # maximal shower length [m]
const double d = sqrt((Z + zmax) * (Z + zmax) + R * R);
d = sqrt((Z + zmax) * (Z + zmax) + R * R);
if (C * t > std::max(n1 * D, zmax + n1 * d)) {
if (C * t > max(n1 * D, zmax + n1 * d)) {
return value;
}
end
if (C * t < n0 * D) {
JRoot rz(R, n0, t + Z / C); // square root
JRoot rz(R, n0, t + Z / C) # square root
if (!rz.is_valid) {
return value;
}
end
if (rz.second > Z) {
zmin = rz.second - Z;
}
end
if (rz.first > Z) {
zmin = rz.first - Z;
}
}
end
end
if (C * t > n1 * D) {
JRoot rz(R, n1, t + Z / C); // square root
JRoot rz(R, n1, t + Z / C) # square root
if (!rz.is_valid) {
return value;
}
end
if (rz.second > Z) {
zmin = rz.second - Z;
}
end
if (rz.first > Z) {
zmin = rz.first - Z;
}
}
end
end
if (C * t < zmax + n0 * d) {
JRoot rz(R, n0, t + Z / C); // square root
JRoot rz(R, n0, t + Z / C) # square root
if (!rz.is_valid) {
return value;
}
end
if (rz.first > Z) {
zmax = rz.first - Z;
}
end
if (rz.second > Z) {
zmax = rz.second - Z;
}
}
end
end
if (zmin < 0.0) {
zmin = 0.0;
}
end
if (zmax > zmin) {
const double ymin = geanz.getIntegral(E, zmin);
const double ymax = geanz.getIntegral(E, zmax);
const double dy = (ymax - ymin) / size();
ymin = geanz.getIntegral(E, zmin);
ymax = geanz.getIntegral(E, zmax);
dy = (ymax - ymin) / size();
if (dy > 2 * std::numeric_limits<double>::epsilon()) {
if dy > 2 * eps()
for (double y = ymin + 0.5 * dy; y < ymax; y += dy) {
const double z = Z + geanz.getLength(E, y);
const double d = sqrt(R * R + z * z);
const double t1 = t + (Z - z) / C - d * getIndexOfRefraction() / C;
z = Z + geanz.getLength(E, y);
d = sqrt(R * R + z * z);
t1 = t + (Z - z) / C - d * getIndexOfRefraction() / C;
value +=
dy * E * getDirectLightFromEMshower(d, -z / d, theta, phi, t1);
}
}
}
end
end
end
return value;
}
end
"""
scatteredlightfromEMshower(params::LMParameters, pmt::PMTModel, E, D, cd, θ, ϕ, Δt)
......@@ -332,73 +332,69 @@ Probability density function for scattered light from EM-shower. Returns
- `Δt`: time difference relative to direct Cherenkov light
"""
function scatteredlightfromEMshower(params::LMParameters, pmt::PMTModel, E, D, cd, θ, ϕ, Δt)
double getScatteredLightFromEMshower(const double E, const double D_m,
const double cd, const double theta,
const double phi,
const double t_ns) const {
double value = 0;
value = 0.0
const double sd = sqrt((1.0 + cd) * (1.0 - cd));
const double D = std::max(D_m, getRmin());
const double R = D * sd; // minimal distance of approach [m]
const double Z = -D * cd;
const double t = D * getIndexOfRefraction() / C + t_ns; // time [ns]
sd = sqrt((1.0 + cd) * (1.0 - cd));
D = max(D_m, getRmin());
R = D * sd # minimal distance of approach [m]
Z = -D * cd;
t = D * getIndexOfRefraction() / C + t_ns # time [ns]
const double n0 = getIndexOfRefractionGroup(wmax);
n0 = getIndexOfRefractionGroup(wmax);
double zmin = 0.0; // minimal shower length [m]
double zmax = geanz.getLength(E, 1.0); // maximal shower length [m]
zmin = 0.0 # minimal shower length [m]
zmax = geanz.getLength(E, 1.0) # maximal shower length [m]
const double d = sqrt((Z + zmax) * (Z + zmax) + R * R);
d = sqrt((Z + zmax) * (Z + zmax) + R * R);
if (C * t < n0 * D) {
JRoot rz(R, n0, t + Z / C); // square root
JRoot rz(R, n0, t + Z / C) # square root
if (!rz.is_valid) {
return value;
}
end
if (rz.second > Z) {
zmin = rz.second - Z;
}
end
if (rz.first > Z) {
zmin = rz.first - Z;
}
}
end
end
if (C * t < zmax + n0 * d) {
JRoot rz(R, n0, t + Z / C); // square root
JRoot rz(R, n0, t + Z / C) # square root
if (!rz.is_valid) {
return value;
}
end
if (rz.first > Z) {
zmax = rz.first - Z;
}
end
if (rz.second > Z) {
zmax = rz.second - Z;
}
}
end
end
const double ymin = geanz.getIntegral(E, zmin);
const double ymax = geanz.getIntegral(E, zmax);
const double dy = (ymax - ymin) / size();
ymin = geanz.getIntegral(E, zmin);
ymax = geanz.getIntegral(E, zmax);
dy = (ymax - ymin) / size();
if (dy > 2 * std::numeric_limits<double>::epsilon()) {
if (dy > 2 * eps()) {
for (double y = ymin + 0.5 * dy; y < ymax; y += dy) {
const double z = Z + geanz.getLength(E, y);
const double d = sqrt(R * R + z * z);
const double t1 = t + (Z - z) / C - d * getIndexOfRefraction() / C;
z = Z + geanz.getLength(E, y);
d = sqrt(R * R + z * z);
t1 = t + (Z - z) / C - d * getIndexOfRefraction() / C;
value +=
dy * E * getScatteredLightFromEMshower(d, -z / d, theta, phi, t1);
}
}
end
end
return value;
}
end
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment