GammaLib 2.0.0
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCOMModelDRBPhibarBins.cpp
Go to the documentation of this file.
1/***************************************************************************
2 * GCOMModelDRBPhibarBins.cpp - COMPTEL DRB model Phibar bin fitting class *
3 * ----------------------------------------------------------------------- *
4 * copyright (C) 2021 by Juergen Knoedlseder *
5 * ----------------------------------------------------------------------- *
6 * *
7 * This program is free software: you can redistribute it and/or modify *
8 * it under the terms of the GNU General Public License as published by *
9 * the Free Software Foundation, either version 3 of the License, or *
10 * (at your option) any later version. *
11 * *
12 * This program is distributed in the hope that it will be useful, *
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15 * GNU General Public License for more details. *
16 * *
17 * You should have received a copy of the GNU General Public License *
18 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19 * *
20 ***************************************************************************/
21/**
22 * @file GCOMModelDRBPhibarBins.cpp
23 * @brief COMPTEL DRB model Phibar bin fitting class implementation
24 * @author Juergen Knoedlseder
25 */
26
27/* __ Includes ___________________________________________________________ */
28#ifdef HAVE_CONFIG_H
29#include <config.h>
30#endif
31#include <typeinfo>
32#include "GException.hpp"
33#include "GTools.hpp"
34#include "GModelRegistry.hpp"
35#include "GEvent.hpp"
36#include "GObservation.hpp"
37#include "GXmlElement.hpp"
39#include "GCOMObservation.hpp"
40#include "GCOMEventBin.hpp"
41
42/* __ Constants __________________________________________________________ */
43
44/* __ Globals ____________________________________________________________ */
46const GModelRegistry g_com_drb_bin_fitting_registry(&g_com_drb_bin_fitting_seed);
47
48/* __ Method name definitions ____________________________________________ */
49#define G_EVAL "GCOMModelDRBPhibarBins::eval(GEvent&, GObservation&, bool&)"
50#define G_NRED "GCOMModelDRBPhibarBins::npred(GEnergy&, GTime&, "\
51 "GObservation&)"
52#define G_MC "GCOMModelDRBPhibarBins::mc(GObservation&, GRan&)"
53#define G_READ "GCOMModelDRBPhibarBins::read(GXmlElement&)"
54#define G_WRITE "GCOMModelDRBPhibarBins::write(GXmlElement&)"
55
56/* __ Macros _____________________________________________________________ */
57
58/* __ Coding definitions _________________________________________________ */
59
60/* __ Debug definitions __________________________________________________ */
61
62
63/*==========================================================================
64 = =
65 = Constructors/destructors =
66 = =
67 ==========================================================================*/
68
69/***********************************************************************//**
70 * @brief Void constructor
71 *
72 * Constructs an empty COMPTEL DRB Phibar bin fitting model.
73 ***************************************************************************/
75{
76 // Initialise members
78
79 // Return
80 return;
81}
82
83
84/***********************************************************************//**
85 * @brief XML constructor
86 *
87 * @param[in] xml XML element.
88 *
89 * Constructs a COMPTEL DRB Phibar bin fitting model from the information
90 * that is found in an XML element. Please refer to the method
91 * GCOMModelDRBPhibarBins::read to learn more about the information that is
92 * expected in the XML element.
93 ***************************************************************************/
95 GModelData(xml)
96{
97 // Initialise members
99
100 // Read XML
101 read(xml);
102
103 // Return
104 return;
105}
106
107
108/***********************************************************************//**
109 * @brief Copy constructor
110 *
111 * @param[in] model COMPTEL DRB Phibar bin fitting model.
112 *
113 * Constructs a COMPTEL DRB Phibar bin fitting model by copying information
114 * from an existing model. Note that the copy is a deep copy, so the original
115 * object can be destroyed after the copy without any loss of information.
116 ***************************************************************************/
118 GModelData(model)
119{
120 // Initialise private members for clean destruction
121 init_members();
122
123 // Copy members
124 copy_members(model);
125
126 // Return
127 return;
128}
129
130
131/***********************************************************************//**
132 * @brief Destructor
133 *
134 * Destroys a COMPTEL DRB Phibar bin fitting model.
135 ***************************************************************************/
137{
138 // Free members
139 free_members();
140
141 // Return
142 return;
143}
144
145
146/*==========================================================================
147 = =
148 = Operators =
149 = =
150 ==========================================================================*/
151
152/***********************************************************************//**
153 * @brief Assignment operator
154 *
155 * @param[in] model COMPTEL DRB Phibar bin fitting model.
156 * @return COMPTEL DRB Phibar bin fitting model
157 *
158 * Assigns the information from a COMPTEL DRB Phibar bin fitting model to the
159 * actual object. Note that a deep copy of the information is performed, so
160 * the original object can be destroyed after the assignment without any loss
161 * of information.
162 ***************************************************************************/
164{
165 // Execute only if object is not identical
166 if (this != &model) {
167
168 // Copy base class members
169 this->GModelData::operator=(model);
170
171 // Free members
172 free_members();
173
174 // Initialise members
175 init_members();
176
177 // Copy members
178 copy_members(model);
179
180 } // endif: object was not identical
181
182 // Return
183 return *this;
184}
185
186
187/*==========================================================================
188 = =
189 = Public methods =
190 = =
191 ==========================================================================*/
192
193/***********************************************************************//**
194 * @brief Clear instance
195 *
196 * Resets the object to a clean initial state. All information that resided
197 * in the object will be lost.
198 ***************************************************************************/
200{
201 // Free class members (base and derived classes, derived class first)
202 free_members();
204 this->GModel::free_members();
205
206 // Initialise members
207 this->GModel::init_members();
209 init_members();
210
211 // Return
212 return;
213}
214
215
216/***********************************************************************//**
217 * @brief Clone instance
218 *
219 * @return Pointer to deep copy of COMPTEL DRB Phibar bin fitting model.
220 *
221 * Clone a COMPTEL DRB Phibar bin fitting model. Cloning performs a deep copy
222 * of the information, so the original object can be destroyed after cloning
223 * without any loss of information.
224 ***************************************************************************/
229
230
231/***********************************************************************//**
232 * @brief Evaluate function
233 *
234 * @param[in] event Observed event.
235 * @param[in] obs Observation.
236 * @param[in] gradients Compute gradients?
237 * @return Background model value.
238 *
239 * @exception GException::invalid_argument
240 * Observation is not a COMPTEL observation.
241 * Event is not a COMPTEL event bin.
242 * @exception GException::invalid_value
243 * Model incompatible with data space.
244 *
245 * Evaluates the COMPTEL DRB Phibar bin fitting model.
246 ***************************************************************************/
248 const GObservation& obs,
249 const bool& gradients) const
250{
251 // Extract COMPTEL observation
252 const GCOMObservation* observation = dynamic_cast<const GCOMObservation*>(&obs);
253 if (observation == NULL) {
254 std::string cls = std::string(typeid(&obs).name());
255 std::string msg = "Observation of type \""+cls+"\" is not a COMPTEL "
256 "observations. Please specify a COMPTEL observation "
257 "as argument.";
259 }
260
261 // Extract COMPTEL event bin
262 const GCOMEventBin* bin = dynamic_cast<const GCOMEventBin*>(&event);
263 if (bin == NULL) {
264 std::string cls = std::string(typeid(&event).name());
265 std::string msg = "Event of type \""+cls+"\" is not a COMPTEL event. "
266 "Please specify a COMPTEL event as argument.";
268 }
269
270 // Get total number of Phibar layers and number of pixels
271 int nphi = observation->drb().map().nmaps();
272 int npix = observation->drb().map().npix();
273 if (nphi != m_values.size()) {
274 std::string msg = "Number of "+gammalib::str(nphi)+" Phibar layers "
275 "differs from number of "+gammalib::str(m_values.size())+
276 " DRB Phibar bin fitting model parameters. Please "
277 "specify a model that is compatible with the data "
278 "space.";
280 }
281 if (npix == 0) {
282 std::string msg = "There are no Chi/Psi pixels in the DRB model. "
283 "Please specify a model that contains Chi/Psi "
284 "pixels.";
286 }
287
288 // Initialise value
289 double value = 0.0;
290
291 // Optionally initialise gradients
292 if (gradients) {
293 for (int i = 0; i < m_values.size(); ++i) {
294 m_values[i].factor_gradient(0.0);
295 }
296 }
297
298 // Get bin size.
299 double size = bin->size();
300
301 // Continue only if bin size is positive
302 if (size > 0.0) {
303
304 // Get bin index and Phibar layer
305 int index = bin->index();
306 int iphi = index / npix;
307
308 // Get DRB model value
309 value = observation->drb().map().pixels()[index] / size;
310
311 // Get scale factor
312 double scale = m_values[iphi].value();
313
314 // Optionally compute gradients
315 if (gradients) {
316
317 // Compute partial derivative
318 double grad = (m_values[iphi].is_free())
319 ? value * m_values[iphi].scale() : 0.0;
320
321 // Set gradient
322 m_values[iphi].factor_gradient(grad);
323
324 } // endif: gradient computation was requested
325
326 // Compute background value
327 value *= scale;
328
329 // Compile option: Check for NaN/Inf
330 #if defined(G_NAN_CHECK)
331 if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
332 std::cout << "*** ERROR: GCOMModelDRBPhibarBins::eval";
333 std::cout << "(index=" << index << "):";
334 std::cout << " NaN/Inf encountered";
335 std::cout << " (iphi=" << iphi;
336 std::cout << ", value=" << value;
337 std::cout << ", scale=" << scale;
338 std::cout << ")" << std::endl;
339 }
340 #endif
341
342 } // endif: binsize was positive
343
344 // Return
345 return value;
346}
347
348
349/***********************************************************************//**
350 * @brief Return spatially integrated data model
351 *
352 * @param[in] obsEng Measured event energy.
353 * @param[in] obsTime Measured event time.
354 * @param[in] obs Observation.
355 * @return Spatially integrated data model.
356 *
357 * @exception GException::feature_not_implemented
358 * Method is not implemented.
359 *
360 * Spatially integrates the data model for a given measured event energy and
361 * event time.
362 *
363 * @todo Implement method.
364 ***************************************************************************/
366 const GTime& obsTime,
367 const GObservation& obs) const
368{
369 // Initialise result
370 double npred = 0.0;
371
372 // Method is not implemented
373 std::string msg = "Spatial integration of data model is not implemented.";
375
376 // Return
377 return npred;
378}
379
380
381/***********************************************************************//**
382 * @brief Return simulated list of events
383 *
384 * @param[in] obs Observation.
385 * @param[in] ran Random number generator.
386 * @return COMPTEL event cube.
387 *
388 * @exception GException::feature_not_implemented
389 * Method is not implemented.
390 *
391 * Draws a sample of events from the COMPTEL DRB Phibar bin fitting model
392 * using a Monte Carlo simulation.
393 *
394 * @todo Implement method.
395 ***************************************************************************/
397 GRan& ran) const
398{
399 // Initialise new event cube
400 GCOMEventCube* cube = new GCOMEventCube;
401
402 // Method is not implemented
403 std::string msg = "Monte Carlo simulation of data model is not implemented.";
405
406 // Return
407 return cube;
408}
409
410
411/***********************************************************************//**
412 * @brief Read model from XML element
413 *
414 * @param[in] xml XML element.
415 *
416 * Read the COMPTEL DRB Phibar bin fitting model from an XML element.
417 *
418 * The model is composed of normalization values that define the scale
419 * factors for all Phibar layers. The following XML file syntax is
420 * expected:
421 *
422 * <source name="Background" type="DRBPhibarBins" instrument="COM">
423 * <parameter name="Normalization" .../>
424 * ...
425 * <parameter name="Normalization" .../>
426 * </source>
427 ***************************************************************************/
429{
430 // Free space for bins
431 m_values.clear();
432
433 // Get number of bins from XML file
434 int bins = xml.elements("parameter");
435
436 // Loop over all bins
437 for (int i = 0; i < bins; ++i) {
438
439 // Allocate normalization parameter
440 GModelPar normalization;
441
442 // Read normalization parameter
443 const GXmlElement* par = xml.element("parameter", i);
444 normalization.read(*par);
445
446 // Set parameter name
447 std::string normalization_name = "Phibar normalization "+gammalib::str(i);
448
449 // Set normalization attributes
450 normalization.name(normalization_name);
451 normalization.unit("");
452 normalization.has_grad(true);
453
454 // Push normalization parameter on list
455 m_values.push_back(normalization);
456
457 } // endfor: looped over bins
458
459 // Read model attributes
460 read_attributes(xml);
461
462 // Set pointers
463 set_pointers();
464
465 // Return
466 return;
467}
468
469
470/***********************************************************************//**
471 * @brief Write model into XML element
472 *
473 * @param[in] xml XML element.
474 *
475 * @exception GException::invalid_value
476 * Existing XML element is not of required type
477 * Invalid number of model parameters or bins found in XML element.
478 *
479 * Write the COMPTEL DRB Phibar bin fitting model into an XML element.
480 *
481 * The model is composed of normalization values that define the scale
482 * factors for all Phibar layers. The following XML file syntax is
483 * expected:
484 *
485 * <source name="Background" type="DRBPhibarBins" instrument="COM">
486 * <parameter name="Normalization" .../>
487 * ...
488 * <parameter name="Normalization" .../>
489 * </source>
490 ***************************************************************************/
492{
493 // Initialise pointer on source
494 GXmlElement* src = NULL;
495
496 // Search corresponding source
497 int n = xml.elements("source");
498 for (int k = 0; k < n; ++k) {
499 GXmlElement* element = xml.element("source", k);
500 if (element->attribute("name") == name()) {
501 src = element;
502 break;
503 }
504 }
505
506 // If no source with corresponding name was found then append one.
507 // Set also the type and the instrument.
508 if (src == NULL) {
509 src = xml.append("source");
510 src->attribute("name", name());
511 src->attribute("type", type());
512 if (instruments().length() > 0) {
513 src->attribute("instrument", instruments());
514 }
515 }
516
517 // Verify model type
518 if (src->attribute("type") != type()) {
519 std::string msg = "Invalid model type \""+src->attribute("type")+"\" "
520 "found in XML file. Model type \""+type()+"\" "
521 "expected.";
523 }
524
525 // Determine number of bins
526 int bins = m_values.size();
527
528 // If XML element has 0 parameters then append parameters
529 if (src->elements() == 0) {
530 for (int i = 0; i < bins; ++i) {
531 src->append(GXmlElement("parameter"));
532 }
533 }
534
535 // Verify that XML element has the required number of bins
536 if (src->elements() != bins || src->elements("parameter") != bins) {
537 std::string msg = "Invalid number of bins "+
538 gammalib::str(src->elements("parameter"))+
539 " found in XML file, but model requires exactly "+
540 gammalib::str(bins)+" bins.";
542 }
543
544 // Loop over all bins
545 for (int i = 0; i < bins; ++i) {
546
547 // Get bin
548 GXmlElement* par = src->element("parameter", i);
549 GModelPar mpar = m_values[i];
550 mpar.name("Normalization");
551 mpar.write(*par);
552
553 } // endfor: looped over bins
554
555 // Write model attributes
556 write_attributes(*src);
557
558 // Return
559 return;
560}
561
562
563/***********************************************************************//**
564 * @brief Print model information
565 *
566 * @param[in] chatter Chattiness.
567 * @return String containing model information.
568 ***************************************************************************/
569std::string GCOMModelDRBPhibarBins::print(const GChatter& chatter) const
570{
571 // Initialise result string
572 std::string result;
573
574 // Continue only if chatter is not silent
575 if (chatter != SILENT) {
576
577 // Append header
578 result.append("=== GCOMModelDRBPhibarBins ===");
579
580 // Append attributes
581 result.append("\n"+print_attributes());
582
583 // Append bins summary
584 result.append("\n"+gammalib::parformat("Number of parameters"));
585 result.append(gammalib::str(size()));
586 for (int i = 0; i < size(); ++i) {
587 result.append("\n"+m_pars[i]->print(chatter));
588 }
589
590 } // endif: chatter was not silent
591
592 // Return result
593 return result;
594}
595
596
597/*==========================================================================
598 = =
599 = Private methods =
600 = =
601 ==========================================================================*/
602
603/***********************************************************************//**
604 * @brief Initialise class members
605 ***************************************************************************/
607{
608 // Initialise members
609 m_values.clear();
610
611 // Return
612 return;
613}
614
615
616/***********************************************************************//**
617 * @brief Copy class members
618 *
619 * @param[in] model Model.
620 ***************************************************************************/
622{
623 // Copy members
624 m_values = model.m_values;
625
626 // Set parameter pointers
627 set_pointers();
628
629 // Return
630 return;
631}
632
633
634/***********************************************************************//**
635 * @brief Delete class members
636 ***************************************************************************/
638{
639 // Return
640 return;
641}
642
643
644/***********************************************************************//**
645 * @brief Set pointers
646 *
647 * Set pointers to all model parameters. The pointers are stored in a vector
648 * that is member of the GModel base class.
649 ***************************************************************************/
651{
652 // Clear parameter pointer(s)
653 m_pars.clear();
654
655 // Get number of bins
656 int bins = m_values.size();
657
658 // Set parameter pointers for all bins
659 for (int i = 0; i < bins; ++i) {
660
661 // Signal that values have gradients
662 m_values[i].has_grad(true);
663
664 // Set pointer
665 m_pars.push_back(&(m_values[i]));
666
667 } // endfor: looped over bins
668
669 // Return
670 return;
671}
#define G_WRITE
COMPTEL event bin class interface definition.
const GCOMModelDRBPhibarBins g_com_drb_bin_fitting_seed
#define G_NRED
COMPTEL DRB model Phibar bin fitting class interface definition.
COMPTEL observation class interface definition.
Abstract event base class definition.
Exception handler interface definition.
#define G_EVAL
Model registry class definition.
Abstract observation base class interface definition.
Gammalib tools definition.
GChatter
Definition GTypemaps.hpp:33
@ SILENT
Definition GTypemaps.hpp:34
XML element node class interface definition.
const GSkyMap & map(void) const
Return DRI sky map.
Definition GCOMDri.hpp:267
COMPTEL event bin class.
const int & index(void) const
Return bin index.
virtual double size(void) const
Return size of event bin.
COMPTEL event bin container class.
COMPTEL DRB model Phibar bin fitting class.
virtual void write(GXmlElement &xml) const
Write model into XML element.
virtual ~GCOMModelDRBPhibarBins(void)
Destructor.
virtual void read(const GXmlElement &xml)
Read model from XML element.
virtual GCOMModelDRBPhibarBins & operator=(const GCOMModelDRBPhibarBins &model)
Assignment operator.
virtual std::string type(void) const
Return model type.
void free_members(void)
Delete class members.
virtual void clear(void)
Clear instance.
void set_pointers(void)
Set pointers.
std::vector< GModelPar > m_values
Normalisation values.
void copy_members(const GCOMModelDRBPhibarBins &model)
Copy class members.
GCOMModelDRBPhibarBins(void)
Void constructor.
virtual GCOMEventCube * mc(const GObservation &obs, GRan &ran) const
Return simulated list of events.
virtual double npred(const GEnergy &obsEng, const GTime &obsTime, const GObservation &obs) const
Return spatially integrated data model.
virtual std::string print(const GChatter &chatter=NORMAL) const
Print model information.
virtual double eval(const GEvent &event, const GObservation &obs, const bool &gradients=false) const
Evaluate function.
virtual GCOMModelDRBPhibarBins * clone(void) const
Clone instance.
void init_members(void)
Initialise class members.
Interface class for COMPTEL observations.
const GCOMDri & drb(void) const
Return background model.
Class that handles energies in a unit independent way.
Definition GEnergy.hpp:48
Abstract interface for the event classes.
Definition GEvent.hpp:71
Abstract data model class.
void init_members(void)
Initialise class members.
void free_members(void)
Delete class members.
virtual GModelData & operator=(const GModelData &model)
Assignment operator.
Model parameter class.
Definition GModelPar.hpp:87
void write(GXmlElement &xml) const
Set or update parameter attributes in XML element.
void read(const GXmlElement &xml)
Extract parameter attributes from XML element.
Interface definition for the model registry class.
GModelPar & scale(const int &index)
Returns reference to scale parameter by index.
Definition GModel.cpp:369
std::vector< GModelPar * > m_pars
Pointers to all model parameters.
Definition GModel.hpp:178
void write_attributes(GXmlElement &xml) const
Write model attributes.
Definition GModel.cpp:723
void free_members(void)
Delete class members.
Definition GModel.cpp:672
void init_members(void)
Initialise class members.
Definition GModel.cpp:622
std::string print_attributes(void) const
Print model attributes.
Definition GModel.cpp:770
int size(void) const
Return number of parameters in model.
Definition GModel.hpp:233
std::string instruments(void) const
Returns instruments to which model applies.
Definition GModel.cpp:310
void read_attributes(const GXmlElement &xml)
Read model attributes.
Definition GModel.cpp:684
const std::string & name(void) const
Return parameter name.
Definition GModel.hpp:261
Abstract observation base class.
bool has_grad(void) const
Signal if parameter gradient is computed analytically.
const std::string & unit(void) const
Return parameter unit.
const std::string & name(void) const
Return parameter name.
Random number generator class.
Definition GRan.hpp:44
const int & nmaps(void) const
Returns number of maps.
Definition GSkyMap.hpp:389
const int & npix(void) const
Returns number of pixels.
Definition GSkyMap.hpp:345
const double * pixels(void) const
Returns pointer to pixel data.
Definition GSkyMap.hpp:475
Time class.
Definition GTime.hpp:55
XML element node class.
const GXmlAttribute * attribute(const int &index) const
Return attribute.
virtual GXmlNode * append(const GXmlNode &node)
Append XML child node.
Definition GXmlNode.cpp:287
virtual GXmlElement * element(const int &index)
Return pointer to GXMLElement child.
Definition GXmlNode.cpp:640
virtual int elements(void) const
Return number of GXMLElement children of node.
Definition GXmlNode.cpp:586
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition GTools.cpp:1143
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition GTools.hpp:184
bool is_notanumber(const double &x)
Signal if argument is not a number.
Definition GTools.hpp:201
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition GTools.cpp:489