Simon Bonnefoy / gammalib

Commit 4104319bea27dc3c08b07e46672e23a1b74018a6

Authored by Jürgen Knödlseder
1 parent 3f8fb819

Complete implementation of general CTA background model (#2694) 

The general CTA background model now correctly works for the existing radial acceptance models. Hence we can start to implement the new model types.
Inline Side-by-side
Showing 4 changed files with 101 additions and 184 deletions
inst/cta/include/GCTAModelRadial.hpp
  View file @4104319
... ... @@ -56,19 +56,15 @@ public:
56 56  
57 57 // Operators
58 58 virtual GCTAModelRadial& operator=(const GCTAModelRadial& model);
59   - virtual GModelPar& operator[](const int& index);
60   - virtual const GModelPar& operator[](const int& index) const;
61   - virtual GModelPar& operator[](const std::string& name);
62   - virtual const GModelPar& operator[](const std::string& name) const;
63 59  
64 60 // Implemented virtual methods
65   - virtual double eval(const GCTAInstDir& dir,
66   - const GEnergy& energy,
67   - const GTime& time,
68   - const bool& gradients = false) const;
69   - virtual GCTAInstDir mc(const GEnergy& energy,
70   - const GTime& time,
71   - GRan& ran) const;
  61 + virtual double eval(const GCTAInstDir& dir,
  62 + const GEnergy& energy,
  63 + const GTime& time,
  64 + const bool& gradients = false) const;
  65 + virtual GCTAInstDir mc(const GEnergy& energy,
  66 + const GTime& time,
  67 + GRan& ran) const;
72 68  
73 69 // Pure virtual methods
74 70 virtual void clear(void) = 0;
... ... @@ -83,31 +79,13 @@ public:
83 79 // Derived class pure virtual methods
84 80 virtual double eval(const double& offset,
85 81 const bool& gradients = false) const = 0;
86   - virtual GCTAInstDir mc(const GCTAInstDir& dir, GRan& ran) const = 0;
87   -
88   - // Methods
89   - int size(void) const;
  82 + virtual GCTAInstDir mc(GRan& ran) const = 0;
90 83  
91 84 protected:
92 85 // Protected methods
93 86 void init_members(void);
94 87 void copy_members(const GCTAModelRadial& model);
95 88 void free_members(void);
96   -
97   - // Proteced members
98   - std::vector<GModelPar*> m_pars; //!< Parameter pointers
99 89 };
100 90  
101   -
102   -/***********************************************************************//**
103   - * @brief Return number of model parameters
104   - *
105   - * @return Number of model parameters.
106   - ***************************************************************************/
107   -inline
108   -int GCTAModelRadial::size(void) const
109   -{
110   - return ((int)m_pars.size());
111   -}
112   -
113 91 #endif /* GCTAMODELRADIAL_HPP */
... ...
inst/cta/pyext/GCTAModelRadial.i
  View file @4104319
... ... @@ -42,20 +42,29 @@ public:
42 42 GCTAModelRadial(const GCTAModelRadial& model);
43 43 virtual ~GCTAModelRadial(void);
44 44  
  45 + // Implemented virtual methods
  46 + virtual double eval(const GCTAInstDir& dir,
  47 + const GEnergy& energy,
  48 + const GTime& time,
  49 + const bool& gradients = false) const;
  50 + virtual GCTAInstDir mc(const GEnergy& energy,
  51 + const GTime& time,
  52 + GRan& ran) const;
  53 +
  54 +
45 55 // Pure virtual methods
46 56 virtual void clear(void) = 0;
47 57 virtual GCTAModelRadial* clone(void) const = 0;
48 58 virtual std::string classname(void) const = 0;
49 59 virtual std::string type(void) const = 0;
50   - virtual double eval(const double& offset,
51   - const bool& gradients = false) const = 0;
52   - virtual GCTAInstDir mc(const GCTAInstDir& dir, GRan& ran) const = 0;
53 60 virtual double omega(void) const = 0;
54 61 virtual void read(const GXmlElement& xml) = 0;
55 62 virtual void write(GXmlElement& xml) const = 0;
56 63  
57   - // Methods
58   - int size(void) const;
  64 + // Derived class pure virtual methods
  65 + virtual double eval(const double& offset,
  66 + const bool& gradients = false) const = 0;
  67 + virtual GCTAInstDir mc(GRan& ran) const = 0;
59 68 };
60 69  
61 70  
... ...
inst/cta/src/GCTAModelBackground.cpp
  View file @4104319
... ... @@ -551,18 +551,26 @@ GCTAEventList* GCTAModelBackground::mc(const GObservation&amp; obs,
551 551 // TODO: Handle special case of cube spatial model. This will
552 552 // replace the spectral model by a spectral nodes model.
553 553  
  554 + // Get solid angle of spatial model. This only works for an energy
  555 + // independent spatial model!!!
  556 + double solidangle = m_spatial->npred(GEnergy(), GTime(), obs);
  557 +
554 558 // Loop over all energy boundaries
555 559 for (int ieng = 0; ieng < ebounds.size(); ++ieng) {
556 560  
557 561 // Compute the background rate in model within the energy
558   - // boundaries from spectral component (units: cts/s).
559   - // Note that the time here is ontime. Deadtime correction will
560   - // be done later.
561   - double rate = spectral->flux(ebounds.emin(ieng), ebounds.emax(ieng));
  562 + // boundaries from spectral component (units: cts/s/sr).
  563 + double flux = spectral->flux(ebounds.emin(ieng), ebounds.emax(ieng));
  564 +
  565 + // Derive expecting rate (units: cts/s). Note that the time here
  566 + // is good time. Deadtime correction will be done later.
  567 + double rate = flux * solidangle;
562 568  
563 569 // Debug option: dump rate
564 570 #if defined(G_DUMP_MC)
565 571 std::cout << "GCTAModelBackground::mc(\"" << name() << "\": ";
  572 + std::cout << "flux=" << flux << " cts/s/sr, ";
  573 + std::cout << "solidangle=" << solidangle << " sr, ";
566 574 std::cout << "rate=" << rate << " cts/s)" << std::endl;
567 575 #endif
568 576  
... ... @@ -593,8 +601,8 @@ GCTAEventList* GCTAModelBackground::mc(const GObservation&amp; obs,
593 601 #if defined(G_DUMP_MC)
594 602 std::cout << " Interval " << itime;
595 603 std::cout << " events=" << n_events << std::endl;
596   - int n_killed_by_deadtime = 0;
597   - int n_killed_by_roi = 0;
  604 + int n_removed_by_deadtime = 0;
  605 + int n_trial_outside_roi = 0;
598 606 #endif
599 607  
600 608 // Loop over events
... ... @@ -605,7 +613,7 @@ GCTAEventList* GCTAModelBackground::mc(const GObservation&amp; obs,
605 613 if (deadc < 1.0) {
606 614 if (ran.uniform() > deadc) {
607 615 #if defined(G_DUMP_MC)
608   - n_killed_by_deadtime++;
  616 + n_removed_by_deadtime++;
609 617 #endif
610 618 continue;
611 619 }
... ... @@ -617,42 +625,53 @@ GCTAEventList* GCTAModelBackground::mc(const GObservation&amp; obs,
617 625 times[i],
618 626 ran);
619 627  
620   - // Get Monte Carlo event direction from spatial model.
621   - // This only will set the DETX and DETY coordinates.
622   - GCTAInstDir instdir = m_spatial->mc(energy, times[i], ran);
  628 + // Get an instrument direction within the RoI. This is
  629 + // potentially tried 100 times so that if we really can't
  630 + // a valid instrument direction the code is not locked up
  631 + for (int k = 0; k < 100; ++k) {
623 632  
624   - // Derive sky direction from instrument coordinates
625   - GSkyDir skydir = pnt.skydir(instdir);
  633 + // Get Monte Carlo event direction from spatial model.
  634 + // This only will set the DETX and DETY coordinates.
  635 + GCTAInstDir instdir = m_spatial->mc(energy, times[i], ran);
626 636  
627   - // Set sky direction in GCTAInstDir object
628   - instdir.dir(skydir);
  637 + // Derive sky direction from instrument coordinates
  638 + GSkyDir skydir = pnt.skydir(instdir);
629 639  
630   - // Allocate event
631   - GCTAEventAtom event;
  640 + // Set sky direction in GCTAInstDir object
  641 + instdir.dir(skydir);
632 642  
633   - // Set event attributes
634   - event.dir(instdir);
635   - event.energy(energy);
636   - event.time(times[i]);
  643 + // Allocate event
  644 + GCTAEventAtom event;
637 645  
638   - // Append event to list if it falls in ROI
639   - if (events.roi().contains(event)) {
640   - list->append(event);
641   - }
642   - #if defined(G_DUMP_MC)
643   - else {
644   - n_killed_by_roi++;
645   - }
646   - #endif
  646 + // Set event attributes
  647 + event.dir(instdir);
  648 + event.energy(energy);
  649 + event.time(times[i]);
  650 +
  651 + // Append event to list if it falls in RoI and break
  652 + // the look
  653 + if (events.roi().contains(event)) {
  654 + list->append(event);
  655 + break;
  656 + }
  657 +
  658 + // ... otherwise optionally bookkeep the trial
  659 + #if defined(G_DUMP_MC)
  660 + else {
  661 + n_trial_outside_roi++;
  662 + }
  663 + #endif
  664 +
  665 + } // endfor: trial look for instrument direction within RoI
647 666  
648 667 } // endfor: looped over all events
649 668  
650   - // Debug option: provide statisics
  669 + // Debug option: provide statisics
651 670 #if defined(G_DUMP_MC)
652   - std::cout << " Killed by deadtime=";
653   - std::cout << n_killed_by_deadtime << std::endl;
654   - std::cout << " Killed by ROI=";
655   - std::cout << n_killed_by_roi << std::endl;
  671 + std::cout << " Events removed due to deadtime=";
  672 + std::cout << n_removed_by_deadtime << std::endl;
  673 + std::cout << " Trials outside RoI=";
  674 + std::cout << n_trial_outside_roi << std::endl;
656 675 #endif
657 676  
658 677 } // endfor: looped over all GTIs
... ... @@ -987,7 +1006,7 @@ void GCTAModelBackground::set_pointers(void)
987 1006 * @param[in] spatial XML element.
988 1007 * @return Pointer to spatial model.
989 1008 *
990   - * Returns pointer to a radial model that is defined in an XML element.
  1009 + * Returns pointer to a spatial model that is defined in an XML element.
991 1010 ***************************************************************************/
992 1011 GCTAModelSpatial* GCTAModelBackground::xml_spatial(const GXmlElement& spatial) const
993 1012 {
... ...
inst/cta/src/GCTAModelRadial.cpp
  View file @4104319
... ... @@ -20,7 +20,7 @@
20 20 ***************************************************************************/
21 21 /**
22 22 * @file GCTAModelRadial.cpp
23   - * @brief Abstract radial acceptance model class implementation
  23 + * @brief Abstract radial background model class implementation
24 24 * @author Juergen Knoedlseder
25 25 */
26 26  
... ... @@ -66,7 +66,7 @@ GCTAModelRadial::GCTAModelRadial(void)
66 66 /***********************************************************************//**
67 67 * @brief Copy constructor
68 68 *
69   - * @param[in] model Radial acceptance model.
  69 + * @param[in] model Radial background model.
70 70 ***************************************************************************/
71 71 GCTAModelRadial::GCTAModelRadial(const GCTAModelRadial& model)
72 72 {
... ... @@ -103,9 +103,9 @@ GCTAModelRadial::~GCTAModelRadial(void)
103 103 /***********************************************************************//**
104 104 * @brief Assignment operator
105 105 *
106   - * @param[in] model Radial acceptance model.
  106 + * @param[in] model Radial background model.
107 107 ***************************************************************************/
108   -GCTAModelRadial& GCTAModelRadial::operator= (const GCTAModelRadial& model)
  108 +GCTAModelRadial& GCTAModelRadial::operator=(const GCTAModelRadial& model)
109 109 {
110 110 // Execute only if object is not identical
111 111 if (this != &model) {
... ... @@ -126,100 +126,6 @@ GCTAModelRadial&amp; GCTAModelRadial::operator= (const GCTAModelRadial&amp; model)
126 126 }
127 127  
128 128  
129   -/***********************************************************************//**
130   - * @brief Returns model parameter
131   - *
132   - * @param[in] index Parameter index [0,...,size()-1].
133   - *
134   - * @exception GException::out_of_range
135   - * Parameter index is out of range.
136   - ***************************************************************************/
137   -GModelPar& GCTAModelRadial::operator[](const int& index)
138   -{
139   - // Compile option: raise exception if index is out of range
140   - #if defined(G_RANGE_CHECK)
141   - if (index < 0 || index >= size())
142   - throw GException::out_of_range(G_ACCESS1, index, 0, size()-1);
143   - #endif
144   -
145   - // Return reference
146   - return *(m_pars[index]);
147   -}
148   -
149   -
150   -/***********************************************************************//**
151   - * @brief Returns model parameter (const version)
152   - *
153   - * @param[in] index Parameter index [0,...,size()-1].
154   - *
155   - * @exception GException::out_of_range
156   - * Parameter index is out of range.
157   - ***************************************************************************/
158   -const GModelPar& GCTAModelRadial::operator[](const int& index) const
159   -{
160   - // Compile option: raise exception if index is out of range
161   - #if defined(G_RANGE_CHECK)
162   - if (index < 0 || index >= size())
163   - throw GException::out_of_range(G_ACCESS1, index, 0, size()-1);
164   - #endif
165   -
166   - // Return reference
167   - return *(m_pars[index]);
168   -}
169   -
170   -
171   -/***********************************************************************//**
172   - * @brief Returns reference to model parameter
173   - *
174   - * @param[in] name Parameter name.
175   - *
176   - * @exception GException::par_not_found
177   - * Parameter with specified name not found in container.
178   - ***************************************************************************/
179   -GModelPar& GCTAModelRadial::operator[](const std::string& name)
180   -{
181   - // Get parameter index
182   - int index = 0;
183   - for (; index < size(); ++index) {
184   - if (m_pars[index]->name() == name)
185   - break;
186   - }
187   -
188   - // Throw exception if parameter name was not found
189   - if (index >= size())
190   - throw GException::par_not_found(G_ACCESS2, name);
191   -
192   - // Return reference
193   - return *(m_pars[index]);
194   -}
195   -
196   -
197   -/***********************************************************************//**
198   - * @brief Returns reference to model parameter (const version)
199   - *
200   - * @param[in] name Parameter name.
201   - *
202   - * @exception GException::par_not_found
203   - * Parameter with specified name not found in container.
204   - ***************************************************************************/
205   -const GModelPar& GCTAModelRadial::operator[](const std::string& name) const
206   -{
207   - // Get parameter index
208   - int index = 0;
209   - for (; index < size(); ++index) {
210   - if (m_pars[index]->name() == name)
211   - break;
212   - }
213   -
214   - // Throw exception if parameter name was not found
215   - if (index >= size())
216   - throw GException::par_not_found(G_ACCESS2, name);
217   -
218   - // Return reference
219   - return *(m_pars[index]);
220   -}
221   -
222   -
223 129 /*==========================================================================
224 130 = =
225 131 = Public methods =
... ... @@ -230,39 +136,50 @@ const GModelPar&amp; GCTAModelRadial::operator[](const std::string&amp; name) const
230 136 * @brief Evaluate function
231 137 *
232 138 * @param[in] dir Event direction.
233   - * @param[in] energy Event energy.
234   - * @param[in] time Event time.
  139 + * @param[in] energy Event energy (not used).
  140 + * @param[in] time Event time (not used).
235 141 * @param[in] gradients Compute gradients?
236 142 * @return Function value
237 143 *
238   - * @todo Implement method
  144 + * Evaluate radial model for a given event direction. The energy and time of
  145 + * the event are not used.
239 146 ***************************************************************************/
240 147 double GCTAModelRadial::eval(const GCTAInstDir& dir,
241 148 const GEnergy& energy,
242 149 const GTime& time,
243 150 const bool& gradients) const
244 151 {
245   - // Return
246   - return 0.0;
  152 + // Compute offset angle in degrees
  153 + double offset = dir.theta() * gammalib::rad2deg;
  154 +
  155 + // Evaluate function
  156 + double value = eval(offset, gradients);
  157 +
  158 + // Return value
  159 + return value;
247 160 }
248 161  
249 162  
250 163 /***********************************************************************//**
251 164 * @brief Returns MC instrument direction
252 165 *
253   - * @param[in] energy Event energy.
254   - * @param[in] time Event time.
  166 + * @param[in] energy Event energy (not used).
  167 + * @param[in] time Event time (not used).
255 168 * @param[in,out] ran Random number generator.
256 169 * @return Instrument direction
257 170 *
258   - * @todo Implement method
  171 + * Return random instrument direction. The energy and time of the event are
  172 + * not used.
259 173 ***************************************************************************/
260 174 GCTAInstDir GCTAModelRadial::mc(const GEnergy& energy,
261 175 const GTime& time,
262 176 GRan& ran) const
263 177 {
264   - // Return
265   - return GCTAInstDir();
  178 + // Get random instrument direction
  179 + GCTAInstDir dir = mc(ran);
  180 +
  181 + // Return instrument direction
  182 + return dir;
266 183 }
267 184  
268 185  
... ... @@ -277,9 +194,6 @@ GCTAInstDir GCTAModelRadial::mc(const GEnergy&amp; energy,
277 194 ***************************************************************************/
278 195 void GCTAModelRadial::init_members(void)
279 196 {
280   - // Initialise members
281   - m_pars.clear();
282   -
283 197 // Return
284 198 return;
285 199 }
... ... @@ -292,9 +206,6 @@ void GCTAModelRadial::init_members(void)
292 206 ***************************************************************************/
293 207 void GCTAModelRadial::copy_members(const GCTAModelRadial& model)
294 208 {
295   - // Copy members
296   - m_pars = model.m_pars;
297   -
298 209 // Return
299 210 return;
300 211 }
... ...