/* 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 .
*
* Copyright Chris Rogers, 2016
*/
#include
#include
#include
#include
#include
#include "gsl/gsl_odeiv.h"
#include "gsl/gsl_errno.h"
#include "Geant4/G4Navigator.hh"
#include "Geant4/G4TransportationManager.hh"
#include "Geant4/G4VPhysicalVolume.hh"
#include "Geant4/G4LogicalVolume.hh"
#include "Geant4/G4Material.hh"
#include "src/common_cpp/Utils/Globals.hh"
#include "src/common_cpp/Utils/Squeak.hh"
#include "src/common_cpp/Simulation/GeometryNavigator.hh"
#include "src/common_cpp/Recon/Global/MaterialModelDynamic.hh"
#include "src/common_cpp/Recon/Global/MaterialModelAxialLookup.hh"
#include "src/common_cpp/Recon/Kalman/Global/ErrorTrackingControl.hh"
#include "src/common_cpp/Recon/Kalman/Global/ErrorTrackingControlLookup.hh"
#include "src/common_cpp/Recon/Kalman/Global/CentroidTracking.hh"
namespace MAUS {
namespace Kalman {
namespace Global {
std::stringstream CentroidTracking::tracking_fail;
const double CentroidTracking::_float_tolerance = 1e-15;
CentroidTracking* CentroidTracking::_tz_for_propagate = NULL;
CentroidTracking::CentroidTracking()
: _field(Globals::GetMCFieldConstructor()), _charge(1.) {
}
CentroidTracking::TrackingSuccess CentroidTracking::Propagate(double x[8], double target_z) {
double test = x[0]+x[1]+x[2]+x[3];
if (test != test) {
Squeak::mout(Squeak::debug) <<
"Position Not A Number in CentroidTracking::Propagate" << std::endl;
return failure;
}
_mass_squared = x[4]*x[4] - x[5]*x[5] - x[6]*x[6] - x[7]*x[7];
if (_mass_squared != _mass_squared) {
Squeak::mout(Squeak::debug) << "Momentum Not a Number in "
<< "CentroidTracking::Propagate" << std::endl;
return failure;
}
if (_mass_squared < -_float_tolerance) {
Squeak::mout(Squeak::debug) << "Negative mass in CentroidTracking::Propagate" << std::endl;
return failure;
}
const gsl_odeiv_step_type * T = gsl_odeiv_step_rk4;
gsl_odeiv_step * step = gsl_odeiv_step_alloc(T, 8);
gsl_odeiv_control * control = NULL;
if (_track_model == em_forwards_dynamic ||
_track_model == em_backwards_dynamic) {
if (_geometry == geant4) {
control = ErrorTrackingControl::gsl_odeiv_control_et_new(
_min_step_size, _max_step_size);
} else if (_geometry == axial_lookup) {
control = ErrorTrackingControlLookup::gsl_odeiv_control_lookup_et_new(
_min_step_size, _max_step_size);
}
}
gsl_odeiv_evolve * evolve = gsl_odeiv_evolve_alloc(8);
gsl_odeiv_system system =
{CentroidTracking::EquationsOfMotion, NULL, 8, NULL};
double z = x[3];
_gsl_h = _min_step_size;
if (z > target_z) {
_gsl_h *= -1; // stepping backwards...
_charge = -1;
} else {
_charge = 1; // stepping forwards...
}
int nsteps = 0;
_tz_for_propagate = this;
tracking_fail.str(""); // clear the error stream
CentroidTracking::TrackingSuccess ret_value = success;
while (fabs(z-target_z) > _float_tolerance) {
nsteps++;
int status = gsl_odeiv_evolve_apply
(evolve, NULL, step, &system, &z, target_z, &_gsl_h, x);
if (control != NULL) {
control->type->hadjust(control->state, 0, 0, x, NULL, NULL, &_gsl_h);
}
if (nsteps > _max_n_steps) {
tracking_fail << "Exceeded maximum number of steps\n";
status = failure;
}
if (status != GSL_SUCCESS) {
tracking_fail << "Killing tracking with step size " <<_gsl_h
<< " at step " << nsteps << " of " << _max_n_steps << "\n";
print(tracking_fail, x);
}
if (status == failure) {
Squeak::mout(Squeak::debug) << tracking_fail.str() << std::endl;
ret_value = failure; // and free memory
break;
}
}
gsl_odeiv_evolve_free(evolve);
gsl_odeiv_control_free(control);
gsl_odeiv_step_free(step);
return ret_value;
}
int CentroidTracking::EquationsOfMotion(double z,
const double x[8],
double dxdz[8],
void* params) {
if (fabs(x[7]) < _float_tolerance) {
// z-momentum is 0
tracking_fail << "pz " << x[7] << " was less than float tolerance "
<< _float_tolerance << "\n";
return GSL_FAILURE;
}
for (size_t i = 0; i < 8; ++i) {
if (x[i] != x[i]) {
tracking_fail << "Detected nan in input i: " << i
<< " value: " << x[i] << "\n";
print(tracking_fail, x);
return GSL_FAILURE;
}
}
if (_tz_for_propagate == NULL) {
tracking_fail << "_tz_for_propagate was NULL\n";
return GSL_FAILURE;
}
int em_success = EMEquationOfMotion(z, x, dxdz, params);
if (em_success != GSL_SUCCESS) {
tracking_fail << "EM transport failed\n";
return em_success;
}
if (_tz_for_propagate->_eloss_model == no_eloss) {
return GSL_SUCCESS;
}
int material_success = MaterialEquationOfMotion(z, x, dxdz, params);
if (material_success != GSL_SUCCESS) {
tracking_fail << "Material transport failed\n";
return material_success;
}
return GSL_SUCCESS;
}
int CentroidTracking::EMEquationOfMotion(double z,
const double x[8],
double dxdz[8],
void* params) {
double field[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
double xfield[4] = {x[1], x[2], x[3], x[0]};
_tz_for_propagate->_field->GetFieldValue(xfield, field);
double direction = -1.;
if ((_tz_for_propagate->_track_model == em_forwards_dynamic ||
_tz_for_propagate->_track_model == em_forwards_static) &&
_tz_for_propagate->_gsl_h < 0.) {
// forwards track propagating in the backwards direction
direction = -1.;
} else if (_tz_for_propagate->_gsl_h > 0.) {
// backwards track propagating in the forwards direction
direction = +1.;
}
double dir = fabs(x[7])/x[7];
double charge = _tz_for_propagate->_charge;
double dtdz = x[4]/x[7];
dxdz[0] = dtdz/c_l; // dt/dz
dxdz[1] = x[5]/x[7]; // dx/dz = px/pz
dxdz[2] = x[6]/x[7]; // dy/dz = py/pz
dxdz[3] = 1.0; // dz/dz
// dE/dz only contains electric field as B conserves energy, not relevant at
// least in step 4 as all fields are static.
dxdz[4] = (dxdz[1]*charge*field[3] + dxdz[2]*charge*field[4] +
charge*field[5])*dir; // dE/dz
// dpx/dz = q*c*(dy/dz*Bz - dz/dz*By) + q*Ex*dt/dz
dxdz[5] = direction*charge*c_l*(dxdz[2]*field[2] - dxdz[3]*field[1])
+ charge*field[3]*dtdz*dir; // dpx/dz
dxdz[6] = direction*charge*c_l*(dxdz[3]*field[0] - dxdz[1]*field[2])
+ charge*field[4]*dtdz*dir; // dpy/dz
dxdz[7] = direction*charge*c_l*(dxdz[1]*field[1] - dxdz[2]*field[0])
+ charge*field[5]*dtdz*dir; // dpz/dz
return GSL_SUCCESS;
}
int CentroidTracking::MaterialEquationOfMotion(double z,
const double x[8],
double dxdz[8],
void* params) {
// Must be called after EMEquationOfMotion
MaterialModel* material = NULL;
if (_tz_for_propagate->_geometry == geant4) {
material = new MaterialModelDynamic(x[1], x[2], x[3]);
} else if (_tz_for_propagate->_geometry == axial_lookup) {
material = new MaterialModelAxialLookup(x[1], x[2], x[3]);
} else {
Squeak::mout(Squeak::warning) << "Did not recognise material model";
return failure;
}
if (material->GetMaterial() == NULL) { // material is not enabled
delete material;
return GSL_SUCCESS;
}
double energy = sqrt(x[4]*x[4]);
double p = sqrt(x[5]*x[5]+x[6]*x[6]+x[7]*x[7]);
double mass = sqrt(energy*energy-p*p);
double charge = _tz_for_propagate->_charge;
double dEdz = 0.;
if (_tz_for_propagate->_eloss_model == most_prob_forwards ||
_tz_for_propagate->_eloss_model == most_prob_backwards) {
material->SetWillUseMeanEnergyLoss(false);
} else {
material->SetWillUseMeanEnergyLoss(true);
}
if (_tz_for_propagate->_eloss_model == bethe_bloch_forwards ||
_tz_for_propagate->_eloss_model == most_prob_forwards) {
dEdz = material->dEdx(energy, mass, charge); // should be negative
} else if (_tz_for_propagate->_eloss_model == bethe_bloch_backwards ||
_tz_for_propagate->_eloss_model == most_prob_backwards) {
// should be negative; during tracking backwards dz is negative so dEdz
// gets handled properly
dEdz = -material->dEdx(energy, mass, charge);
} else {
dEdz = 0.; // no eloss
}
// E^2 = p^2+m^2; 2E dE/dz = 2p dp/dz
double dpdz = energy/p*dEdz;
dxdz[4] += dEdz;
dxdz[5] += dpdz*x[5]/p;
dxdz[6] += dpdz*x[6]/p;
dxdz[7] += dpdz*x[7]/p;
delete material;
return GSL_SUCCESS;
}
void CentroidTracking::SetField(const BTField* field) {
if (field == NULL) {
_field = Globals::GetMCFieldConstructor();
} else {
_field = field;
}
}
void CentroidTracking::SetGeometryModel(std::string geometry_model) {
if (geometry_model == "geant4") {
SetGeometryModel(geant4);
} else if (geometry_model == "axial_lookup") {
SetGeometryModel(axial_lookup);
} else {
throw Exceptions::Exception(Exceptions::recoverable,
"Failed to recognise geometry model '"+geometry_model
+"' Should be 'geant4' or 'axial_lookup'",
"CentroidTracking::SetGeometryModel");
}
}
std::ostream& CentroidTracking::print(std::ostream& out, const double* x) {
out << "t: " << x[0] << " pos: " << x[1] << " " << x[2] << " " << x[3] << "\n"
<< "E: " << x[4] << " mom: " << x[5] << " " << x[6] << " " << x[7] << std::endl;
return out;
}
void CentroidTracking::SetMinStepSize(double min_step_size) {
// step size is always positive (direction is determined automatically)
min_step_size = fabs(min_step_size);
if (min_step_size > _max_step_size) {
throw Exceptions::Exception(Exceptions::recoverable,
"Min step size should be < max step size",
"CentroidTracking::SetMinStepSize");
}
_min_step_size = min_step_size;
}
void CentroidTracking::SetMaxStepSize(double max_step_size) {
// step size is always positive (direction is determined automatically)
max_step_size = fabs(max_step_size);
if (_min_step_size > max_step_size) {
throw Exceptions::Exception(Exceptions::recoverable,
"Min step size should be < max step size",
"CentroidTracking::SetMaxStepSize");
}
_max_step_size = max_step_size;
}
void CentroidTracking::SetTrackingModel(std::string tracking_model) {
if (tracking_model == "em_forwards_dynamic") {
SetTrackingModel(em_forwards_dynamic);
} else if (tracking_model == "em_backwards_dynamic") {
SetTrackingModel(em_backwards_dynamic);
} else if (tracking_model == "em_forwards_static") {
SetTrackingModel(em_forwards_static);
} else if (tracking_model == "em_backwards_static") {
SetTrackingModel(em_backwards_static);
} else {
throw Exceptions::Exception(
Exceptions::recoverable,
"Did not recognise tracking model "+tracking_model,
"CentroidTracking::SetTrackingModel");
}
}
} // namespace Global
} // namespace Kalman
} // namespace MAUS