Iftach Sadeh / ctools

Commit 0c21591f85a335f55d954691c76d83fbb864881c

Authored by Luigi Tibaldo
1 parent b0f2c752

running version, does nothing, methods not tested yet

Inline Side-by-side
Showing 3 changed files with 56 additions and 4 deletions
cscripts/Makefile.am
  View file @0c21591
... ... @@ -33,6 +33,7 @@ pythonscripts = $(top_builddir)/cscripts/__init__.py \
33 33 $(srcdir)/csobsinfo.py \
34 34 $(srcdir)/csobsselect.py \
35 35 $(srcdir)/csphasecrv.py \
  36 + $(srcdir)/csphagen.py \
36 37 $(srcdir)/cspull.py \
37 38 $(srcdir)/csresmap.py \
38 39 $(srcdir)/csroot2caldb.py \
... ... @@ -64,6 +65,7 @@ parfiles = $(srcdir)/cscaldb.par \
64 65 $(srcdir)/csobsdef.par \
65 66 $(srcdir)/csobsinfo.par \
66 67 $(srcdir)/csobsselect.par \
  68 + $(srcdir)/csphagen.par \
67 69 $(srcdir)/csphasecrv.par \
68 70 $(srcdir)/cspull.par \
69 71 $(srcdir)/csresmap.par \
... ... @@ -93,6 +95,7 @@ linkscripts = cscaldb \
93 95 csobsdef \
94 96 csobsinfo \
95 97 csobsselect \
  98 + csphagen \
96 99 csphasecrv \
97 100 cspull \
98 101 csresmap \
... ...
cscripts/csphagen.par
  View file @0c21591
... ... @@ -37,7 +37,7 @@ dec, r, a, 22.01,-90,90, "Declination of source region centre (deg)"
37 37 glon, r, a, 184.56,0,360, "Galactic longitude of source region centre (deg)"
38 38 glat, r, a, -5.79,-90,90, "Galactic latitude of source region centre (deg)"
39 39 rad, r, a, 0.2,0,180, "Radius of source region circle (deg)"
40   -bkgregmin, i, h, 2,1,, "Minimum number of background region for REFLECTED method"
  40 +bkgregmin, i, h, 2,1,1000, "Minimum number of background region for REFLECTED method"
41 41 stack, b, a, yes,,, "Stack all runs into single pha (if no produce run-by-run files)"
42 42  
43 43 #
... ...
cscripts/csphagen.py
  View file @0c21591
... ... @@ -21,6 +21,24 @@
21 21 # ==========================================================================
22 22 import gammalib
23 23 import ctools
  24 +import math
  25 +
  26 +
  27 +def atan2(y, x):
  28 + """
  29 + helper function to solve triangulation problem
  30 + """
  31 + if x > 0:
  32 + val = math.atan(y / x)
  33 + elif x < 0 and y >= 0:
  34 + val = math.atan(y / x) + math.pi
  35 + elif x < 0 and y < 0:
  36 + val = math.atan(y / x) - math.pi
  37 + elif x == 0 and y > 0:
  38 + val = math.pi / 2
  39 + elif x == 0 and y < 0:
  40 + val = math.pi / 2
  41 + return val
24 42  
25 43  
26 44 # =============== #
... ... @@ -37,7 +55,7 @@ class csphagen(ctools.cscript):
37 55 Constructor
38 56 """
39 57 # Set name, version
40   - self._name = 'csfindobs'
  58 + self._name = 'csphagen'
41 59 self._version = ctools.__version__
42 60  
43 61 # Initialise observation container from constructor arguments
... ... @@ -90,12 +108,17 @@ class csphagen(ctools.cscript):
90 108 gammalib.GSkyRegionCircle(self._src_dir, self._rad))
91 109  
92 110 # Stacking
93   - self._stack = self['stack'].bool()
  111 + self._stack = self['stack'].boolean()
94 112  
95 113 # Query ahead output parameters
96 114 if (self._read_ahead()):
97 115 self['outroot'].string()
98 116  
  117 + # If there are no observations in container then get them from the
  118 + # parameter file
  119 + if self._obs.size() == 0:
  120 + self._obs = self._get_observations(False)
  121 +
99 122 # Write input parameters into logger
100 123 self._log_parameters(gammalib.TERSE)
101 124  
... ... @@ -112,6 +135,32 @@ class csphagen(ctools.cscript):
112 135 offset = pnt_dir.dist_deg(self._src_dir)
113 136 outregions = []
114 137  
  138 + if self._srcshape == "CIRCLE":
  139 + # angular separation of reflected regions wrt camera center
  140 + # and number
  141 + alpha = 2 * math.asin(2 * self._rad / offset)
  142 + N = int(2 * math.pi) / alpha
  143 + alpha = 2 * math.pi / N
  144 + # loop to create reflected regions
  145 + for s in range(2, N - 1):
  146 + phi = s * alpha
  147 + # region center
  148 + xii = offset(1 - 2 * math.sin(phi / 2))
  149 + yii = 2 * offset * math.sqrt(1 - math.sin(phi / 2) ** 2)
  150 + if phi > math.pi / 2:
  151 + yii *= 1
  152 + theta = atan2(pnt_dir.dec_deg() - self._src_dir.dec_deg(),
  153 + pnt_dir.ra_deg() - self._src_dir.ra_deg())
  154 + xi = xii * math.cos(theta) - yii * math.sin(theta)
  155 + yi = xii * math.sin(theta) + yii * math.cos(theta)
  156 + x = xi + pnt_dir.ra_deg()
  157 + y = yi + pnt_dir.dec_deg()
  158 + ctr_dir = gammalib.GSkyDir()
  159 + ctr_dir.radec_deg(x, y)
  160 + # region creation
  161 + region = gammalib.GSkyRegionCircle(ctr_dir, self._rad)
  162 + outregions.append(region)
  163 +
115 164 return outregions
116 165  
117 166 # Public methods
... ... @@ -140,7 +189,7 @@ class csphagen(ctools.cscript):
140 189 for obs in self._obs:
141 190 bkg_reg = gammalib.GSkyRegions()
142 191 if self._bkgmethod == "REFLECTED":
143   - regions = self.reflected_regions(obs)
  192 + regions = self._reflected_regions(obs)
144 193 for region in regions:
145 194 bkg_reg.append(region)
146 195 onoff = gammalib.GCTAOnOffObservation(obs, etrue, ereco,
... ...