/* This file is part of MAUS: http://micewww.pp.rl.ac.uk/projects/maus
*
* MAUS is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* MAUS is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with MAUS. If not, see .
*
*/
#include
#include
#include "Geant4/G4Material.hh"
#include "src/common_cpp/Simulation/GeometryNavigator.hh"
#include "src/common_cpp/Utils/Squeak.hh"
#include "src/common_cpp/Recon/Global/MaterialModel.hh"
namespace MAUS {
std::map MaterialModel::_enabled_materials;
bool MaterialModel::_default_enablement = true;
void MaterialModel::SetMaterial(const G4Material* material) {
_material = material;
_n_e = _material->GetElectronDensity();
_I = _material->GetIonisation()->GetMeanExcitationEnergy();
_x_0 = _material->GetIonisation()->GetX0density();
_x_1 = _material->GetIonisation()->GetX1density();
_C = _material->GetIonisation()->GetCdensity();
_a = _material->GetIonisation()->GetAdensity();
_k = _material->GetIonisation()->GetMdensity();
_rad_len_ratio = 13.6*13.6/_material->GetRadlen();
_density = _material->GetDensity()/g*cm3;
}
double MaterialModel::MeandEdx(double E, double m, double charge) {
// bethe bloch formula, see PDG particle data book
// assume density effect is 0
double beta2 = (E*E-m*m)/E/E; // beta^2 = p^2/E^2
double gamma = E/m;
double bg2 = beta2*gamma*gamma;
double mRatio = _m_e/m;
double T_max = 2.0*_m_e*bg2/(1.0 + 2.0*gamma*mRatio + mRatio*mRatio);
double coefficient = std::log(2.0*_m_e*bg2*T_max/(_I*_I))/2. - beta2;
// 0.307075 is cm^2 per mol; _density is g/cm^3
coefficient *= 0.307075*charge*charge/beta2*_density/cm;
double dEdx = 0.;
for (size_t i = 0; i < _material->GetNumberOfElements(); ++i) {
double frac = _material->GetFractionVector()[i];
double a_el = _material->GetElement(i)->GetA()/(g/mole);
double z_el = _material->GetElement(i)->GetZ();
dEdx += frac*z_el/a_el*coefficient;
}
return -dEdx;
}
double MaterialModel::MostProbabledEdx(double E, double m, double charge) {
double dens_thickness = _thickness*_density/cm; // g/cm^2
double beta2 = (E*E-m*m)/E/E; // beta^2 = p^2/E^2
double gamma = E/m;
double bg2 = beta2*gamma*gamma;
double meanZOverA = 0;
for (size_t i = 0; i < _material->GetNumberOfElements(); ++i) {
meanZOverA += _material->GetFractionVector()[i]*
_material->GetElement(i)->GetZ()/
(_material->GetElement(i)->GetA()/(g/mole));
}
double xi = 0.307075/2.*meanZOverA*dens_thickness/beta2;
double dE = -xi*(log(2*m*bg2/_I)+log(xi/_I)-beta2); // ignore density effect
return dE/_thickness;
}
double MaterialModel::dEdx(double E, double m, double charge) {
if (_material == NULL) {
return 0.;
}
if (_will_use_mean_energy_loss) {
return MeandEdx(E, m, charge);
} else {
return MostProbabledEdx(E, m, charge);
}
}
double MaterialModel::d2EdxdE(double E, double m, double charge) {
double delta_energy_pos = dEdx(E+_d2EdxdE_delta_const, m, charge);
double delta_energy_neg = dEdx(E-_d2EdxdE_delta_const, m, charge);
double deriv = (delta_energy_pos-delta_energy_neg)/2./_d2EdxdE_delta_const;
return deriv;
}
double MaterialModel::dtheta2dx(double E, double m, double charge) {
// moliere scattering formula, see PDG particle data book
// assume log term is 0
double one_over_beta_c_p = E/(E*E-m*m); // E/p^2
double theta_0_sq = _rad_len_ratio*one_over_beta_c_p*one_over_beta_c_p;
return theta_0_sq;
}
double MaterialModel::estrag2(double E, double m, double charge) {
// NIM A 532, Neuffer, "Introduction to Ionisation Cooling", 2004
// Scraped from Ann. Rev. Nucl. Sci. 13, Fano, 1963
// Note that I ignore the units in the NIM paper as I think it is
// wrong; I think the correct units are MeV^2/cm, not MeV^2 cm^2/g
// (Rogers, 2016)
double gamma2 = E*E/m/m;
double beta2 = (E*E-m*m)/E/E;
double estrag2 = _estrag_const*gamma2*(1-beta2/2.)/cm; // missing factor Z/A
return estrag2;
}
bool MaterialModel::IsEnabledMaterial(std::string material_name) {
bool found = _enabled_materials.find(material_name) !=
_enabled_materials.end();
if (found) {
return _enabled_materials[material_name];
} else if (material_name.size() > 2 && material_name.substr(0, 3) == std::string("G4_")) {
material_name = material_name.substr(3);
}
if (_enabled_materials.find(material_name) != _enabled_materials.end()) {
return _enabled_materials[material_name];
}
return _default_enablement;
}
void MaterialModel::EnableMaterial(std::string material_name) {
_enabled_materials[material_name] = true;
}
void MaterialModel::ClearEnabledMaterials() {
_enabled_materials = std::map();
}
void MaterialModel::SetDefaultMaterialEnablement(bool default_enablement) {
_default_enablement = default_enablement;
}
void MaterialModel::DisableMaterial(std::string material_name) {
_enabled_materials[material_name] = false;
}
void MaterialModel::SetWillUseMeanEnergyLoss(bool will_use_mean_energy_loss) {
_will_use_mean_energy_loss = will_use_mean_energy_loss;
}
bool MaterialModel::GetWillUseMeanEnergyLoss() const {
return _will_use_mean_energy_loss;
}
std::ostream& MaterialModel::PrintMaterials(std::ostream& out) {
out << "Materials:\n";
for (std::map::iterator it = _enabled_materials.begin();
it != _enabled_materials.end(); ++it) {
out << " " << it->first << " " << it->second << "\n";
}
return out;
}
} // namespace MAUS