Interface for the COMPTEL instrument response representation class. More...

#include <GCOMIaq.hpp>

Inheritance diagram for GCOMIaq:

Classes
class	response_kernel

class	smearing_kernel

Public Member Functions
	GCOMIaq (void)
	Void constructor. More...

	GCOMIaq (const GCOMIaq &iaq)
	Copy constructor. More...

	GCOMIaq (const double &phigeo_max, const double &phigeo_bin_size, const double &phibar_max, const double &phibar_bin_size)
	COMPTEL instrument response representation constructor. More...

	~GCOMIaq (void)
	Destructor. More...

GCOMIaq &	operator= (const GCOMIaq &iaq)
	Assignment operator. More...

void	clear (void)
	Clear instance. More...

GCOMIaq *	clone (void) const
	Clone instance. More...

std::string	classname (void) const
	Return class name. More...

void	set (const GEnergy &energy, const GEbounds &ebounds)
	Set mono-energetic IAQ. More...

void	set (const GModelSpectral &spectrum, const GEbounds &ebounds)
	Set continuum IAQ. More...

void	save (const GFilename &filename, const bool &clobber=false) const
	Save COMPTEL instrument response representation into FITS file. More...

std::string	print (const GChatter &chatter=NORMAL) const
	Print COMPTEL instrument response representation information. More...

Public Member Functions inherited from GBase
virtual	~GBase (void)
	Destructor. More...

Private Member Functions
void	init_members (void)
	Initialise class members. More...

void	copy_members (const GCOMIaq &iaq)
	Copy class members. More...

void	free_members (void)
	Delete class members. More...

void	init_response (void)
	Initialise COMPTEL response function and instrument characteristics. More...

void	remove_cards (void)
	Remove any extra header cards. More...

void	compton_kinematics (const double &energy)
	Generate IAQ matrix based on Compton kinematics. More...

double	compute_iaq_bin (const double &etrue1, const double &etrue2, const double &phibar)
	Computes the IAQ for one bin. More...

void	klein_nishina (const double &energy)
	Multiply IAQ matrix by the Klein-Nishina formula. More...

double	klein_nishina_bin (const double &energy, const double &phigeo)
	Computes Klein-Nishina probability for one bin. More...

double	klein_nishina_integral (const double &v, const double &a)
	Computes Klein-Nishina integral. More...

void	weight_iaq (const double &energy)
	Weight IAQ matrix. More...

void	location_smearing (const double &zenith)
	Perform location smearing. More...

std::vector< double >	location_spread (const double &zenith) const
	Compute location spread vector. More...

Private Attributes
GFitsImageFloat	m_iaq
	Response. More...

GEbounds	m_ebounds
	Energy boundaries. More...

GCOMD1Response	m_response_d1
	D1 module response. More...

GCOMD2Response	m_response_d2
	D2 module response. More...

GCOMInstChars	m_ict
	Instrument characteristics. More...

double	m_phigeo_max
	Maximum geometrical scatter angle (deg) More...

double	m_phibar_max
	Maximum Compton scatter angle (deg) More...

double	m_phigeo_bin_size
	Bin size in geometrical scatter angle (deg) More...

double	m_phibar_bin_size
	Bin size in Compton scatter angle (deg) More...

double	m_phibar_resolution
	Bin size for oversampling (deg) More...

double	m_e1min
	Minimum D1 energy (MeV) More...

double	m_e1max
	Maximum D1 energy (MeV) More...

double	m_e2min
	Minimum D2 energy (MeV) More...

double	m_e2max
	Maximum D2 energy (MeV) More...

int	m_num_energies
	Number of energies for continuum IAQ. More...

bool	m_psd_correct
	PSD correction usage flag. More...

double	m_zenith
	Zenith angle for location smearing (deg) More...

Detailed Description

Interface for the COMPTEL instrument response representation class.

Definition at line 52 of file GCOMIaq.hpp.

Constructor & Destructor Documentation

GCOMIaq::GCOMIaq ( void )

Void constructor.

Creates an empty COMPTEL instrument response representation.

Definition at line 79 of file GCOMIaq.cpp.

References init_members().

Referenced by clone().

GCOMIaq::GCOMIaq ( const GCOMIaq & iaq )

Copy constructor.

Parameters

[in] iaq COMPTEL instrument response representation.

Definition at line 94 of file GCOMIaq.cpp.

References copy_members(), and init_members().

GCOMIaq::GCOMIaq	(	const double &	phigeo_max,
		const double &	phigeo_bin_size,
		const double &	phibar_max,
		const double &	phibar_bin_size
	)

COMPTEL instrument response representation constructor.

Parameters

[in]	phigeo_max	Maximum geometrical scatter angle (deg).
[in]	phigeo_bin_size	Bin size in geometrical scatter angle (deg).
[in]	phibar_max	Maximum Compton scatter angle (deg).
[in]	phibar_bin_size	Bin size in Compton scatter angle (deg).

Todo:: Test input argument validity.

Definition at line 117 of file GCOMIaq.cpp.

References GFitsHDU::card(), init_members(), m_iaq, m_phibar_bin_size, m_phibar_max, m_phigeo_bin_size, and m_phigeo_max.

GCOMIaq::~GCOMIaq ( void )

Destructor.

Destroys instance of COMPTEL response object.

Definition at line 162 of file GCOMIaq.cpp.

References free_members().

Member Function Documentation

std::string GCOMIaq::classname ( void ) const

inlinevirtual

Return class name.

Returns: String containing the class name ("GCOMIaq").

Implements GBase.

Definition at line 186 of file GCOMIaq.hpp.

void GCOMIaq::clear ( void )

virtual

Clear instance.

Clears COMPTEL instrument response representation by resetting all members to an initial state. Any information that was present in the object before will be lost.

Implements GBase.

Definition at line 218 of file GCOMIaq.cpp.

References free_members(), and init_members().

GCOMIaq * GCOMIaq::clone ( void ) const

virtual

Clone instance.

Returns: Pointer to deep copy of COMPTEL instrument response representation.

Implements GBase.

Definition at line 237 of file GCOMIaq.cpp.

References GCOMIaq().

void GCOMIaq::compton_kinematics ( const double & energy )

private

Generate IAQ matrix based on Compton kinematics.

Parameters

[in] energy Input photon energy (MeV).

Generates an IAQ matrix based on Compton kinematics only. See the compute_iaq_bin() method for the formula used to compute each bin of the IAQ matrix.

The Phibar dimension of the IAQ is sampled at a resolution that is set by the m_phibar_resolution member.

The code implemented is based on the COMPASS RESPSIT2 function IAQWEI.F (release 1.0, 24-FEB-93).

Definition at line 686 of file GCOMIaq.cpp.

References gammalib::com_energy2(), compute_iaq_bin(), gammalib::deg2rad, m_iaq, m_phibar_bin_size, m_phibar_resolution, m_phigeo_bin_size, GFitsImage::naxes(), and sum().

Referenced by set().

double GCOMIaq::compute_iaq_bin	(	const double &	etrue1,
		const double &	etrue2,
		const double &	phibar
	)

private

Computes the IAQ for one bin.

Parameters

[in]	etrue1	True D1 energy deposit (MeV).
[in]	etrue2	True D2 energy deposit (MeV).
[in]	phibar	Compton scatter angle (deg).

Returns: Response for one IAQ bin

Compute the integral

\[ R(\bar{\varphi}|\hat{E_1},\hat{E_2}) = \int_{E_1^{\rm min}}^{E_1^{\rm max}} R(E_1,E_2|\hat{E_1},\hat{E_2}) dE_1 \]

where \(R(E_1,E_2|\hat{E_1},\hat{E_2})\) is computed using GCOMIaq::response_kernel::eval,

\(E_1\) and \(E_2\) are the measured D1 and D2 energy deposits, respectively, \(\hat{E_1}\) and \(\hat{E_2}\) are the true D1 and D2 energy deposits, respectively, and \(\bar{\varphi}\) the Compton scatter angle.

The code implementation is based on the COMPASS RESPSIT2 functions RESPSC.F (release 2.0, 14-Mar-91) and RESINT.F (release 3.0, 21-Oct-91). Since RESPSC.F only defines the threshold while RESINT.F does the real job, both functions were merged.

Using a continuum response for 3-10 MeV the following range for the number of integration kernel calls and integration iterations were determined:

   eps = 1.0e-4    calls = 17-8139      iter = 5-14
   eps = 1.0e-5    calls = 17-65537     iter = 5-17
   eps = 1.0e-6    calls = 17-524289    iter = 5-20
   eps = 1.0e-7    calls = 17-524289    iter = 5-21
   eps = 1.0e-8    calls = 17-524289    iter = 5-21

Graphically there was no difference between the results for any of the values, so in order to have safety margin, a value of 1.0e-6 was selected. Note that RESINT.F implemented a trapezoid rule with 100 steps.

Definition at line 811 of file GCOMIaq.cpp.

References GCOMD1Response::emax(), GEbounds::emax(), GCOMD1Response::emin(), GEbounds::emin(), GCOMD1Response::ewidth(), m_e1max, m_e1min, m_e2max, m_e2min, m_ebounds, m_response_d1, m_response_d2, GEnergy::MeV(), GCOMD1Response::position(), and GCOMD1Response::sigma().

Referenced by compton_kinematics().

void GCOMIaq::copy_members ( const GCOMIaq & iaq )

private

Copy class members.

Parameters

[in] iaq COMPTEL instrument response representation.

Definition at line 585 of file GCOMIaq.cpp.

References m_e1max, m_e1min, m_e2max, m_e2min, m_ebounds, m_iaq, m_ict, m_num_energies, m_phibar_bin_size, m_phibar_max, m_phibar_resolution, m_phigeo_bin_size, m_phigeo_max, m_psd_correct, m_response_d1, m_response_d2, and m_zenith.

Referenced by GCOMIaq(), and operator=().

void GCOMIaq::free_members ( void )

private

Delete class members.

Definition at line 616 of file GCOMIaq.cpp.

Referenced by clear(), operator=(), and ~GCOMIaq().

void GCOMIaq::init_members ( void )

private

Initialise class members.

< Default: 0.25 deg

< Default: 70 keV

< Default: 20 MeV

< Default: 650 keV

< Default: 30 MeV

< Default: 50 input energies

< Default: use PSD correction

< Default: 25 deg

Definition at line 552 of file GCOMIaq.cpp.

References GCOMInstChars::clear(), GCOMD1Response::clear(), GCOMD2Response::clear(), GFitsImageFloat::clear(), GEbounds::clear(), m_e1max, m_e1min, m_e2max, m_e2min, m_ebounds, m_iaq, m_ict, m_num_energies, m_phibar_bin_size, m_phibar_max, m_phibar_resolution, m_phigeo_bin_size, m_phigeo_max, m_psd_correct, m_response_d1, m_response_d2, and m_zenith.

Referenced by clear(), GCOMIaq(), and operator=().

void GCOMIaq::init_response ( void )

private

Initialise COMPTEL response function and instrument characteristics.

Definition at line 626 of file GCOMIaq.cpp.

References m_ict, m_response_d1, and m_response_d2.

Referenced by set().

void GCOMIaq::klein_nishina ( const double & energy )

private

Multiply IAQ matrix by the Klein-Nishina formula.

Parameters

[in] energy Input photon energy (MeV).

The code implemented is based on the COMPASS RESPSIT2 function IAQWEI.F (release 1.0, 24-FEB-93).

Definition at line 886 of file GCOMIaq.cpp.

References klein_nishina_bin(), m_iaq, m_phigeo_bin_size, GFitsImage::naxes(), and sum().

Referenced by set().

double GCOMIaq::klein_nishina_bin	(	const double &	energy,
		const double &	phigeo
	)

private

Computes Klein-Nishina probability for one bin.

Parameters

[in]	energy	Input photon energy (MeV).
[in]	phigeo	Geometric scatter angle (deg).

Returns: Klein-Nishina probability.

Compute the probability that a photon of energy is Compton scattered in a given \(\varphi_{\rm geo}\) bin using

\[ P = \frac{\int_{\varphi_{\rm geo,min}}^{\varphi_{\rm geo,max}} \sigma_{\rm KN}(E, \varphi_{\rm geo}) d\varphi_{\rm geo}} {\int_{0}^{\pi} \sigma_{\rm KN}(E, \varphi_{\rm geo}) d\varphi_{\rm geo}} \]

where \(\sigma_{\rm KN}(E, \varphi_{\rm geo})\) is the Klein-Nishina cross section.

The code implementation is based on the COMPASS RESPSIT2 function RESPSG.F (release 2.0, 26-Apr-90).

Definition at line 955 of file GCOMIaq.cpp.

References cos(), gammalib::deg2rad, klein_nishina_integral(), m_phigeo_bin_size, and gammalib::mec2.

Referenced by klein_nishina().

double GCOMIaq::klein_nishina_integral	(	const double &	v,
		const double &	a
	)

private

Computes Klein-Nishina integral.

Parameters

[in]	v	Integration limit.
[in]	a	Normalized energy.

Returns: Integrated Klein-Nishina cross section.

Computes the indefinite integral of the Klein-Nishina cross section evaluated at the integration limit v.

The code implementation is based on the COMPASS RESPSIT2 function RESSLV.F (release 1.0, 23-Nov-88). The formula used in that code is the following:

W = V*( 1.D0/A + 1.D0 )**2 W = W + ( 2.D0/(A**2) + 2.D0/A - 1.D0 )* DLOG(1.D0 - V) W = W - (V**2)/2.D0 W = W + 1.D0/(A**2) * 1.D0/(1.D0 - V)

It has been check that the original formula gives the same results as the version that was implemented.

Definition at line 1023 of file GCOMIaq.cpp.

References log().

Referenced by klein_nishina_bin().

void GCOMIaq::location_smearing ( const double & zenith )

private

Perform location smearing.

Parameters

[in] zenith Zenith angle of source (deg).

The code implementation is based on the COMPASS RESPSIT2 function LOCSPR.F (release 1.0, 05-JAN-93).

Definition at line 1376 of file GCOMIaq.cpp.

References GNodeArray::append(), GIntegral::eps(), location_spread(), m_iaq, m_phibar_bin_size, GFitsImage::naxes(), GIntegral::romberg(), and GIntegral::silent().

Referenced by set().

std::vector< double > GCOMIaq::location_spread ( const double & zenith ) const

private

Compute location spread vector.

Parameters

[in] zenith Zenith angle of source (deg).

Returns: Vector of location spreads in Gaussian sigma

The method computes the geometrical smearing widths according to D1 and D2 X,Y,Z location uncertainties. The algorithm is based on a constant D1, D2 location uncertainty for all energies. This is also why the method can be called at the end of the IAQ generation: the location spread does not have to be added for each energy separately.

The location spread is a (mild) function of phigeo and source zenith angle.

The code implementation is based on the COMPASS RESPSIT2 function LOCSPR.F (release 1.0, 05-JAN-93).

Definition at line 1487 of file GCOMIaq.cpp.

References cos(), gammalib::deg2rad, m_iaq, m_phigeo_bin_size, GFitsImage::naxes(), gammalib::rad2deg, sin(), and sqrt().

Referenced by location_smearing().

GCOMIaq & GCOMIaq::operator= ( const GCOMIaq & iaq )

Assignment operator.

Parameters

[in] iaq COMPTEL instrument response representation.

Returns: COMPTEL instrument response representation.

Definition at line 184 of file GCOMIaq.cpp.

References copy_members(), free_members(), and init_members().

std::string GCOMIaq::print ( const GChatter & chatter = NORMAL ) const

virtual

Print COMPTEL instrument response representation information.

Parameters

[in] chatter Chattiness.

Returns: String containing COMPTEL instrument response representation information.

Implements GBase.

Definition at line 469 of file GCOMIaq.cpp.

References GEbounds::emax(), GEbounds::emin(), m_e1max, m_e1min, m_e2max, m_e2min, m_ebounds, m_iaq, m_num_energies, m_phibar_bin_size, m_phibar_max, m_phibar_resolution, m_phigeo_bin_size, m_phigeo_max, m_psd_correct, m_zenith, GFitsImage::naxes(), GFitsImage::naxis(), GFitsImage::npix(), gammalib::parformat(), GEnergy::print(), SILENT, gammalib::str(), and sum().

void GCOMIaq::remove_cards ( void )

private

Remove any extra header cards.

Remove all header cards that may have been attached by the set() methods.

Definition at line 644 of file GCOMIaq.cpp.

References GFitsHDU::has_card(), GFitsHDU::header(), and m_iaq.

Referenced by set().

void GCOMIaq::save	(	const GFilename &	filename,
		const bool &	clobber = `false`
	)		const

Save COMPTEL instrument response representation into FITS file.

Parameters

[in]	filename	FITS file name.
[in]	clobber	Overwrite existing file?

Saves the COMPTEL instrument response representation into a FITS file. If the file exists already and the clobber parameter is true, the method will overwrite the content of the existing file. Otherwise, an exception is thrown.

Definition at line 426 of file GCOMIaq.cpp.

References GFits::append(), GFitsHDU::card(), GEbounds::emax(), GEbounds::emin(), GFitsHDU::extname(), GFilename::extname(), GFilename::has_extname(), m_e1max, m_e1min, m_ebounds, m_iaq, m_phibar_resolution, m_psd_correct, m_zenith, GEnergy::MeV(), GFits::saveto(), and GFilename::url().

void GCOMIaq::set	(	const GEnergy &	energy,
		const GEbounds &	ebounds
	)

Set mono-energetic IAQ.

Parameters

[in]	energy	Input photon energy.
[in]	ebounds	Boundaries of observed energy.

The code implemented is based on the COMPASS RESPSIT2 function SPCIAQ.F (release 1.0, 18-DEC-92).

Todo:: Implement geometrical smearing.

Definition at line 254 of file GCOMIaq.cpp.

References GFitsHDU::card(), compton_kinematics(), init_response(), klein_nishina(), location_smearing(), m_ebounds, m_iaq, m_zenith, GEnergy::MeV(), remove_cards(), and weight_iaq().

void GCOMIaq::set	(	const GModelSpectral &	spectrum,
		const GEbounds &	ebounds
	)

Set continuum IAQ.

Parameters

[in]	spectrum	Input spectrum.
[in]	ebounds	Boundaries of observed energy.

Computes the continuum IAQ based on an input spectrum. The method computes the line IAQ for a m_num_energies logarithmically spaced energies within an energy range that conceivably covers the observed energy band.

The code implementation is loosly based on the COMPASS RESPSIT2 functions SPCIAQ.F, GENWEI.F and DERFAN.F (release 1.0, 18-DEC-92). The code was much simplified.