Domenico Tiziani / gammalib

Commit 0201b1fa7dd67db9c4579e5a7ebe0f444de1392b

Authored by knodlseder
1 parent 455ed488

Add deadtime correction to 

GCTAModelRadialAcceptance::eval,
GCTAModelRadialAcceptance::eval_gradients,
GCTAModelRadialAcceptance::npred, and
 GCTAModelRadialAcceptance::mc methods.
Inline Side-by-side
Showing 1 changed file with 125 additions and 32 deletions
inst/cta/src/GCTAModelRadialAcceptance.cpp
  View file @0201b1f
1 1 /***************************************************************************
2 2 * GCTAModelRadialAcceptance.cpp - Radial acceptance model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2012 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -79,6 +79,8 @@ const GModelRegistry g_cta_radial_acceptance_registry(&g_cta_radial_a
79 79  
80 80 /***********************************************************************//**
81 81 * @brief Void constructor
  82 + *
  83 + * Constructs an empty CTA radial acceptance model.
82 84 ***************************************************************************/
83 85 GCTAModelRadialAcceptance::GCTAModelRadialAcceptance(void) : GModelData()
84 86 {
... ... @@ -94,6 +96,11 @@ GCTAModelRadialAcceptance::GCTAModelRadialAcceptance(void) : GModelData()
94 96 * @brief Constructor
95 97 *
96 98 * @param[in] xml XML element.
  99 + *
  100 + * Constructs a CTA radial acceptance model from the information that is
  101 + * found in a XML element. Please refer to the method
  102 + * GCTAModelRadialAcceptance::read
  103 + * to learn more about the information that is expected in the XML element.
97 104 ***************************************************************************/
98 105 GCTAModelRadialAcceptance::GCTAModelRadialAcceptance(const GXmlElement& xml) :
99 106 GModelData(xml)
... ... @@ -117,6 +124,11 @@ GCTAModelRadialAcceptance::GCTAModelRadialAcceptance(const GXmlElement& xml) :
117 124 *
118 125 * @param[in] radial Radial CTA model component.
119 126 * @param[in] spectral Spectral model component.
  127 + *
  128 + * Constructs a CTA radial acceptance model from a radial and a spectral
  129 + * model component. The temporal component is assumed to be constant.
  130 + * Please refer to the classes GCTAModelRadial and GModelSpectral to learn
  131 + * more about the definition of the radial and spectral components.
120 132 ***************************************************************************/
121 133 GCTAModelRadialAcceptance::GCTAModelRadialAcceptance(const GCTAModelRadial& radial,
122 134 const GModelSpectral& spectral)
... ... @@ -145,6 +157,10 @@ GCTAModelRadialAcceptance::GCTAModelRadialAcceptance(const GCTAModelRadial& radi
145 157 * @brief Copy constructor
146 158 *
147 159 * @param[in] model Radial acceptance model.
  160 + *
  161 + * Constructs a CTA radial acceptance model by copying information from an
  162 + * existing model. Note that the copy is a deep copy, so the original object
  163 + * can be destroyed after the copy without any loss of information.
148 164 ***************************************************************************/
149 165 GCTAModelRadialAcceptance::GCTAModelRadialAcceptance(const GCTAModelRadialAcceptance& model) :
150 166 GModelData(model)
... ... @@ -165,6 +181,8 @@ GCTAModelRadialAcceptance::GCTAModelRadialAcceptance(const GCTAModelRadialAccept
165 181  
166 182 /***********************************************************************//**
167 183 * @brief Destructor
  184 + *
  185 + * Destroys a CTA radial acceptance model.
168 186 ***************************************************************************/
169 187 GCTAModelRadialAcceptance::~GCTAModelRadialAcceptance(void)
170 188 {
... ... @@ -186,6 +204,11 @@ GCTAModelRadialAcceptance::~GCTAModelRadialAcceptance(void)
186 204 * @brief Assignment operator
187 205 *
188 206 * @param[in] model Radial acceptance model.
  207 + *
  208 + * Assigns the information from a CTA radial acceptance model to the actual
  209 + * object. Note that a deep copy of the information is performed, so the
  210 + * original object can be destroyed after the assignment without any loss of
  211 + * information.
189 212 ***************************************************************************/
190 213 GCTAModelRadialAcceptance& GCTAModelRadialAcceptance::operator= (const GCTAModelRadialAcceptance& model)
191 214 {
... ... @@ -219,7 +242,10 @@ GCTAModelRadialAcceptance& GCTAModelRadialAcceptance::operator= (const GCTAModel
219 242  
220 243 /***********************************************************************//**
221 244 * @brief Clear instance
222   -***************************************************************************/
  245 + *
  246 + * Resets the object to a clean initial state. All information that resided
  247 + * in the object will be lost.
  248 + ***************************************************************************/
223 249 void GCTAModelRadialAcceptance::clear(void)
224 250 {
225 251 // Free class members (base and derived classes, derived class first)
... ... @@ -239,7 +265,11 @@ void GCTAModelRadialAcceptance::clear(void)
239 265  
240 266 /***********************************************************************//**
241 267 * @brief Clone instance
242   -***************************************************************************/
  268 + *
  269 + * Clone a CTA radial acceptance model. Cloning performs a deep copy of the
  270 + * information, so the original object can be destroyed after cloning without
  271 + * any loss of information.
  272 + ***************************************************************************/
243 273 GCTAModelRadialAcceptance* GCTAModelRadialAcceptance::clone(void) const
244 274 {
245 275 return new GCTAModelRadialAcceptance(*this);
... ... @@ -252,6 +282,11 @@ GCTAModelRadialAcceptance* GCTAModelRadialAcceptance::clone(void) const
252 282 * @param[in] event Observed event.
253 283 * @param[in] obs Observation.
254 284 *
  285 + * Evaluates tha CTA radial acceptance model which is a factorization of a
  286 + * spatial, spectral and temporal model component. This method also applies
  287 + * a deadtime correction factor, so that the normalization of the model is
  288 + * a real rate (counts/exposure time).
  289 + *
255 290 * @todo Add bookkeeping of last value and evaluate only if argument
256 291 * changed
257 292 * @todo Verify that CTA instrument direction pointer is valid, or better,
... ... @@ -285,6 +320,9 @@ double GCTAModelRadialAcceptance::eval(const GEvent& event,
285 320 // Compute value
286 321 double value = rad * spec * temp;
287 322  
  323 + // Apply deadtime correction
  324 + value *= obs.deadc(event.time());
  325 +
288 326 // Return
289 327 return value;
290 328 }
... ... @@ -299,6 +337,12 @@ double GCTAModelRadialAcceptance::eval(const GEvent& event,
299 337 * @exception GCTAException::no_pointing
300 338 * No valid CTA pointing found in observation
301 339 *
  340 + * Evaluates tha CTA radial acceptance model and parameter gradients. The CTA
  341 + * radial acceptance model is a factorization of a spatial, spectral and
  342 + * temporal model component. This method also applies a deadtime correction
  343 + * factor, so that the normalization of the model is a real rate
  344 + * (counts/exposure time).
  345 + *
302 346 * @todo Add bookkeeping of last value and evaluate only if argument
303 347 * changed
304 348 * @todo Verify that CTA instrument direction pointer is valid, or better,
... ... @@ -314,7 +358,7 @@ double GCTAModelRadialAcceptance::eval_gradients(const GEvent& event,
314 358 // Extract CTA pointing direction
315 359 GCTAPointing* pnt = dynamic_cast<GCTAPointing*>(obs.pointing());
316 360 if (pnt == NULL) {
317   - throw GCTAException::no_pointing(G_EVAL);
  361 + throw GCTAException::no_pointing(G_EVAL_GRADIENTS);
318 362 }
319 363  
320 364 // Get instrument direction
... ... @@ -332,6 +376,9 @@ double GCTAModelRadialAcceptance::eval_gradients(const GEvent&amp; event,
332 376 // Compute value
333 377 double value = rad * spec * temp;
334 378  
  379 + // Apply deadtime correction
  380 + value *= obs.deadc(event.time());
  381 +
335 382 // Multiply factors to radial gradients
336 383 if (radial() != NULL) {
337 384 double fact = spec * temp;
... ... @@ -376,8 +423,9 @@ double GCTAModelRadialAcceptance::eval_gradients(const GEvent&amp; event,
376 423 * @exception GCTAException::no_pointing
377 424 * No valid CTA pointing found in observation
378 425 *
379   - * Spatially integrates the data model for a given measured event energy
380   - * and event time.
  426 + * Spatially integrates the data model for a given measured event energy and
  427 + * event time. This method also applies a deadtime correction factor, so that
  428 + * the normalization of the model is a real rate (counts/exposure time).
381 429 ***************************************************************************/
382 430 double GCTAModelRadialAcceptance::npred(const GEnergy& obsEng,
383 431 const GTime& obsTime,
... ... @@ -424,6 +472,9 @@ double GCTAModelRadialAcceptance::npred(const GEnergy&amp; obsEng,
424 472 npred *= spectral()->eval(obsEng);
425 473 npred *= temporal()->eval(obsTime);
426 474  
  475 + // Apply deadtime correction
  476 + npred *= obs.deadc(obsTime);
  477 +
427 478 } // endif: model was valid
428 479  
429 480 // Return
... ... @@ -438,13 +489,18 @@ double GCTAModelRadialAcceptance::npred(const GEnergy&amp; obsEng,
438 489 * @param[in] ran Random number generator.
439 490 *
440 491 * @exception GCTAException::no_pointing
441   - * No pointing found
  492 + * No CTA pointing found in observation.
  493 + *
  494 + * Draws a sample of events from the radial acceptance model using a Monte
  495 + * Carlo simulation. The pointing information, the energy boundaries and the
  496 + * good time interval for the sampling will be extracted from the observation
  497 + * argument that is passed to the method. The method also requires a random
  498 + * number generator of type GRan which is passed by reference, hence the
  499 + * state of the random number generator will be changed by the method.
442 500 *
443   - * This method requires that
444   - * the pointing,
445   - * the energy boundaries, and
446   - * the good time interval
447   - * of the observation have been set up previously.
  501 + * The method also applies a deadtime correction using a Monte Carlo process,
  502 + * taking into account temporal deadtime variations. For this purpose, the
  503 + * method makes use of the time dependent GObservation::deadc method.
448 504 ***************************************************************************/
449 505 GCTAEventList* GCTAModelRadialAcceptance::mc(const GObservation& obs,
450 506 GRan& ran) const
... ... @@ -475,7 +531,8 @@ GCTAEventList* GCTAModelRadialAcceptance::mc(const GObservation&amp; obs,
475 531 // Compute solid angle used for normalization
476 532 double area = radial()->omega();
477 533  
478   - // Derive expecting rate (units: cts/s)
  534 + // Derive expecting rate (units: cts/s). Note that the time here
  535 + // is good time. Deadtime correction will be done later.
479 536 double rate = flux * area;
480 537  
481 538 // Debug option: dump rate
... ... @@ -490,17 +547,26 @@ GCTAEventList* GCTAModelRadialAcceptance::mc(const GObservation&amp; obs,
490 547 for (int itime = 0; itime < obs.events()->gti().size(); ++itime) {
491 548  
492 549 // Get event arrival times from temporal model
493   - GTimes times = m_temporal->mc(rate,
  550 + GTimes times = m_temporal->mc(rate,
494 551 obs.events()->gti().tstart(itime),
495 552 obs.events()->gti().tstop(itime), ran);
496 553  
497 554 // Reserve space for events
498   - if (times.size() > 0)
  555 + if (times.size() > 0) {
499 556 list->reserve(times.size());
  557 + }
500 558  
501 559 // Loop over events
502 560 for (int i = 0; i < times.size(); ++i) {
503 561  
  562 + // Apply deadtime correction
  563 + double deadc = obs.deadc(times[i]);
  564 + if (deadc < 1.0) {
  565 + if (ran.uniform() > deadc) {
  566 + continue;
  567 + }
  568 + }
  569 +
504 570 // Set event direction
505 571 GCTAInstDir dir = radial()->mc(pnt_dir, ran);
506 572  
... ... @@ -589,6 +655,8 @@ void GCTAModelRadialAcceptance::read(const GXmlElement&amp; xml)
589 655 * @brief Write model into XML element
590 656 *
591 657 * @param[in] xml XML element.
  658 + *
  659 + * @todo Document method.
592 660 ***************************************************************************/
593 661 void GCTAModelRadialAcceptance::write(GXmlElement& xml) const
594 662 {
... ... @@ -647,6 +715,8 @@ void GCTAModelRadialAcceptance::write(GXmlElement&amp; xml) const
647 715  
648 716 /***********************************************************************//**
649 717 * @brief Print model information
  718 + *
  719 + * @todo Document method.
650 720 ***************************************************************************/
651 721 std::string GCTAModelRadialAcceptance::print(void) const
652 722 {
... ... @@ -666,40 +736,48 @@ std::string GCTAModelRadialAcceptance::print(void) const
666 736 result.append("\n"+parformat("Instruments"));
667 737 if (m_instruments.size() > 0) {
668 738 for (int i = 0; i < m_instruments.size(); ++i) {
669   - if (i > 0)
  739 + if (i > 0) {
670 740 result.append(", ");
  741 + }
671 742 result.append(m_instruments[i]);
672 743 }
673 744 }
674   - else
  745 + else {
675 746 result.append("all");
  747 + }
676 748  
677 749 // Append model type
678 750 result.append("\n"+parformat("Model type"));
679 751 if (n_radial > 0) {
680 752 result.append("\""+radial()->type()+"\"");
681   - if (n_spectral > 0 || n_temporal > 0)
  753 + if (n_spectral > 0 || n_temporal > 0) {
682 754 result.append(" * ");
  755 + }
683 756 }
684 757 if (n_spectral > 0) {
685 758 result.append("\""+spectral()->type()+"\"");
686   - if (n_temporal > 0)
  759 + if (n_temporal > 0) {
687 760 result.append(" * ");
  761 + }
688 762 }
689   - if (n_temporal > 0)
  763 + if (n_temporal > 0) {
690 764 result.append("\""+temporal()->type()+"\"");
  765 + }
691 766  
692 767 // Append parameters
693 768 result.append("\n"+parformat("Number of parameters")+str(size()));
694 769 result.append("\n"+parformat("Number of radial par's")+str(n_radial));
695   - for (int i = 0; i < n_radial; ++i)
  770 + for (int i = 0; i < n_radial; ++i) {
696 771 result.append("\n"+(*radial())[i].print());
  772 + }
697 773 result.append("\n"+parformat("Number of spectral par's")+str(n_spectral));
698   - for (int i = 0; i < n_spectral; ++i)
  774 + for (int i = 0; i < n_spectral; ++i) {
699 775 result.append("\n"+(*spectral())[i].print());
  776 + }
700 777 result.append("\n"+parformat("Number of temporal par's")+str(n_temporal));
701   - for (int i = 0; i < n_temporal; ++i)
  778 + for (int i = 0; i < n_temporal; ++i) {
702 779 result.append("\n"+(*temporal())[i].print());
  780 + }
703 781  
704 782 // Return result
705 783 return result;
... ... @@ -714,6 +792,8 @@ std::string GCTAModelRadialAcceptance::print(void) const
714 792  
715 793 /***********************************************************************//**
716 794 * @brief Initialise class members
  795 + *
  796 + * @todo Document method.
717 797 ***************************************************************************/
718 798 void GCTAModelRadialAcceptance::init_members(void)
719 799 {
... ... @@ -731,6 +811,8 @@ void GCTAModelRadialAcceptance::init_members(void)
731 811 * @brief Copy class members
732 812 *
733 813 * @param[in] model Model.
  814 + *
  815 + * @todo Document method.
734 816 ***************************************************************************/
735 817 void GCTAModelRadialAcceptance::copy_members(const GCTAModelRadialAcceptance& model)
736 818 {
... ... @@ -749,6 +831,8 @@ void GCTAModelRadialAcceptance::copy_members(const GCTAModelRadialAcceptance&amp; mo
749 831  
750 832 /***********************************************************************//**
751 833 * @brief Delete class members
  834 + *
  835 + * @todo Document method.
752 836 ***************************************************************************/
753 837 void GCTAModelRadialAcceptance::free_members(void)
754 838 {
... ... @@ -788,16 +872,19 @@ void GCTAModelRadialAcceptance::set_pointers(void)
788 872 if (n_pars > 0) {
789 873  
790 874 // Gather radial parameter pointers
791   - for (int i = 0; i < n_radial; ++i)
  875 + for (int i = 0; i < n_radial; ++i) {
792 876 m_pars.push_back(&((*radial())[i]));
  877 + }
793 878  
794 879 // Gather spectral parameters
795   - for (int i = 0; i < n_spectral; ++i)
  880 + for (int i = 0; i < n_spectral; ++i) {
796 881 m_pars.push_back(&((*spectral())[i]));
  882 + }
797 883  
798 884 // Gather temporal parameters
799   - for (int i = 0; i < n_temporal; ++i)
  885 + for (int i = 0; i < n_temporal; ++i) {
800 886 m_pars.push_back(&((*temporal())[i]));
  887 + }
801 888  
802 889 }
803 890  
... ... @@ -844,12 +931,14 @@ GCTAModelRadial* GCTAModelRadialAcceptance::xml_radial(const GXmlElement&amp; radial
844 931 GCTAModelRadial* ptr = registry.alloc(type);
845 932  
846 933 // If model if valid then read model from XML file
847   - if (ptr != NULL)
  934 + if (ptr != NULL) {
848 935 ptr->read(radial);
  936 + }
849 937  
850 938 // ... otherwise throw an exception
851   - else
  939 + else {
852 940 throw GCTAException::model_invalid_radial(G_XML_RADIAL, type);
  941 + }
853 942  
854 943 // Return pointer
855 944 return ptr;
... ... @@ -876,12 +965,14 @@ GModelSpectral* GCTAModelRadialAcceptance::xml_spectral(const GXmlElement&amp; spect
876 965 GModelSpectral* ptr = registry.alloc(type);
877 966  
878 967 // If model if valid then read model from XML file
879   - if (ptr != NULL)
  968 + if (ptr != NULL) {
880 969 ptr->read(spectral);
  970 + }
881 971  
882 972 // ... otherwise throw an exception
883   - else
  973 + else {
884 974 throw GException::model_invalid_spectral(G_XML_SPECTRAL, type);
  975 + }
885 976  
886 977 // Return pointer
887 978 return ptr;
... ... @@ -908,12 +999,14 @@ GModelTemporal* GCTAModelRadialAcceptance::xml_temporal(const GXmlElement&amp; tempo
908 999 GModelTemporal* ptr = registry.alloc(type);
909 1000  
910 1001 // If model if valid then read model from XML file
911   - if (ptr != NULL)
  1002 + if (ptr != NULL) {
912 1003 ptr->read(temporal);
  1004 + }
913 1005  
914 1006 // ... otherwise throw an exception
915   - else
  1007 + else {
916 1008 throw GException::model_invalid_temporal(G_XML_TEMPORAL, type);
  1009 + }
917 1010  
918 1011 // Return pointer
919 1012 return ptr;
... ...