/* 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 "src/legacy/BeamTools/BTField.hh"
#include "src/legacy/BeamTools/BTFieldConstructor.hh"
#include "src/common_cpp/Recon/Kalman/KalmanTrack.hh"
#ifndef recon_kalman_GlobalCentroidTracking_hh_
#define recon_kalman_GlobalCentroidTracking_hh_
class G4Material;
class BTField;
namespace MAUS {
class MaterialModel;
namespace Kalman {
namespace Global {
/** @class CentroidTracking Propagate a track and accompanying error through
* materials and fields
*
* Track propagation is done by standard Lorentz force law F = q (v ^ B).
*
* Through materials, energy loss yields loss in mean energy according to e.g.
* Bethe Bloch model.
*
* We have to be careful to consider the case where forward-going events are
* extrapolated forwards; forward-going tracks are extrapolated backwards; and
* backward-going tracks are extrapolated backwards or forwards. These different
* cases are handled by the ELossModel and TrackingModel.
*
* Numerical integration is done using GSL RK4 routines. Note that these
* routines require a global function pointer. So during tracking we
* assign a static version of CentroidTracking and a static member function
* pointer to point at the tracking.
*
* ELossModel controls model of energy loss through material
* - bethe_bloch_forwards: dE/dz is negative: energy will decrease if step
* size is positive, irrespective of pz direction
* - bethe_bloch_backwards: dE/dz is negative: energy will increase if step
* size is positive, irrespective of pz direction
* - no_eloss: dE/dz is 0.
* TrackingModel: controls how tracking through the electromagnetic fields is
* done
* - em_forwards_dynamic: track forwards, using dynamic step size routine to
* determine step size (where we choose smaller step size near materials)
* - em_backwards_dynamic: track backwards, using dynamic step size routine to
* determine step size (where we choose smaller step size near materials)
* - em_forwards_static: track forwards, using fixed step size
* - em_backwards_static: track forwards, using fixed step size
* GeometryModel: controls how geometry is managed in the fields
* - geant4: use the full geant4 material model. At each step, lookup up the
* distance to the nearest boundary and use this to determine step size
* - axial_lookup: use a lookup table of materials that is determined at run
* time by looking up the materials on axis. Assumes all materials are
* cylindrical.
*
* Other things:-
* - min_step_size, max_step_size: stepping starts at min_step_size.
* - max_n_steps: maximum number of steps. If tracking goes for more steps,
* tracking will abort with GSL_FAILURE.
*
* TODO: hard coded max_n_steps, no "error tolerance" in dynamic step size
* routines
*/
class CentroidTracking {
public:
enum ELossModel {bethe_bloch_forwards, bethe_bloch_backwards,
most_prob_forwards, most_prob_backwards, no_eloss};
enum TrackingModel {em_forwards_dynamic, em_backwards_dynamic,
em_forwards_static, em_backwards_static};
enum GeometryModel {geant4, axial_lookup};
enum TrackingSuccess {success = GSL_SUCCESS, failure = GSL_FAILURE};
/** Constructor; just sets some default things
*
* Defaults field to Globals::GetMCFieldConstructor()
*/
CentroidTracking();
/** Propagate a track to z position target_z
*
* @param x: array of length 29; first 8 parameters are
* (time, x, y, z), (total energy, px, py, pz)
* x is a borrowed pointer; caller owns memory allocated to x.
* Propagation is done "in-place" i.e. x will be updated.
* @param target_z target position to which the track will be propagated
*/
TrackingSuccess Propagate(double x[8], double target_z);
/** Set the field for use by CentroidTracking
* @param BTField* field:- Set the field used for tracking - CentroidTracking
* cannot handle electric fields properly. If NULL, uses MC field
* defined in Globals. This is a borrowed pointer. Caller
* maintains ownership of the memory allocated to field.
*/
void SetField(const BTField* field);
/** Get the field used by CentroidTracking. Note this is a borrowed pointer,
* borrowed from caller to SetField (Globals::MCFieldConstructor by default)
*/
const BTField* GetField() const {return _field;}
/** Set the maximum step size used by the tracking. Will always take the
* absolute value (tracking determines direction).
*
* @throws Exception if _min_step_size > max_step_size
*/
void SetMaxStepSize(double max_step_size);
/** Get the maximum step size used by the tracking */
double GetMaxStepSize() const {return _max_step_size;}
/** Set the minimum step size used by the tracking. Will always take the
* absolute value (tracking determines direction).
*
* @throws Exception if min_step_size > _max_step_size
*/
void SetMinStepSize(double min_step_size);
/** Get the minimum step size used by the tracking */
double GetMinStepSize() const {return _min_step_size;}
/** Set the energy loss model
* @param eloss_model: string from "bethe_bloch_forwards",
* "bethe_bloch_backwards", "most_probable_forwards",
* "most_probable_backwards", "no_eloss"
*/
void SetEnergyLossModel(std::string eloss_model);
/** Set the energy loss model */
void SetEnergyLossModel(ELossModel eloss_model) {_eloss_model = eloss_model;}
/** Get the energy loss model */
ELossModel GetEnergyLossModel() const {return _eloss_model;}
/** Set the tracking model
* @param tracking_model: string from "em_forwards_dynamic",
* "em_backwards_dynamic", "em_forwards_static",
* "em_backwards_static"
*/
void SetTrackingModel(std::string tracking_model);
/** Set the tracking model */
void SetTrackingModel(TrackingModel track_model) {_track_model = track_model;}
/** Get the tracking model */
TrackingModel GetTrackingModel() const {return _track_model;}
/** Set the geometry model
* @param geometry_model: string from "geant4", "axial_lookup"
*/
void SetGeometryModel(std::string geometry_model);
/** Set the geometry model */
void SetGeometryModel(GeometryModel geometry) {_geometry = geometry;}
/** Get the geometry model */
GeometryModel GetGeometryModel() const {return _geometry;}
/** Set the charge of the particle to be tracked, in units of electron charge
*/
void SetCharge(double charge) {_charge = charge;}
/** Get the charge of the particle to be tracked, in units of electron charge
*/
double GetCharge() const {return _charge;}
/** Print an output/input of the Propagate function to out
*
* @param out print to this output stream
* @param x: array formatted as per Propagate function input
*/
static ostream& print(std::ostream& out, const double x[8]);
private:
CentroidTracking(const CentroidTracking& tracking); // disable copy constructor
// x_in is a vector of size 8, going like t,x,y,z; E,px,py,pz
static int EquationsOfMotion(double z, const double x[8], double dxdt[8],
void* params);
static int EMEquationOfMotion(double z, const double x[8], double dxdt[8],
void* params);
static int MaterialEquationOfMotion(double z, const double x[8], double dxdt[8],
void* params);
const BTField* _field;
// scale for derivative calculation
double _min_step_size = 0.1; // magnitude of the step size
double _max_step_size = 100.; // magnitude of the step size
ELossModel _eloss_model = no_eloss;
TrackingModel _track_model = em_forwards_dynamic;
GeometryModel _geometry = geant4;
int _max_n_steps = 100000;
static constexpr double c_l = 299.792458; // mm*ns^{-1}
// transient...
double _gsl_h = 1.; // the step that gsl uses (can be negative)
double _charge;
double _mass_squared;
static CentroidTracking* _tz_for_propagate;
static const double _float_tolerance;
static std::stringstream tracking_fail;
}; // class CentroidTracking
} // namespace Global
} // namespace Kalman
} // namespace MAUS
#endif