1 | /* ======================================================================== *\
|
---|
2 | !
|
---|
3 | ! *
|
---|
4 | ! * This file is part of MARS, the MAGIC Analysis and Reconstruction
|
---|
5 | ! * Software. It is distributed to you in the hope that it can be a useful
|
---|
6 | ! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
|
---|
7 | ! * It is distributed WITHOUT ANY WARRANTY.
|
---|
8 | ! *
|
---|
9 | ! * Permission to use, copy, modify and distribute this software and its
|
---|
10 | ! * documentation for any purpose is hereby granted without fee,
|
---|
11 | ! * provided that the above copyright notice appear in all copies and
|
---|
12 | ! * that both that copyright notice and this permission notice appear
|
---|
13 | ! * in supporting documentation. It is provided "as is" without express
|
---|
14 | ! * or implied warranty.
|
---|
15 | ! *
|
---|
16 | !
|
---|
17 | !
|
---|
18 | ! Author(s): Thomas Bretz, 07/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
|
---|
19 | !
|
---|
20 | ! Copyright: MAGIC Software Development, 2000-2007
|
---|
21 | !
|
---|
22 | !
|
---|
23 | \* ======================================================================== */
|
---|
24 |
|
---|
25 | /////////////////////////////////////////////////////////////////////////////
|
---|
26 | //
|
---|
27 | // MPointingDevCalc
|
---|
28 | //
|
---|
29 | // Calculates the pointing deviation from the starguider information.
|
---|
30 | //
|
---|
31 | // There are some quality parameters to steer which are the valid
|
---|
32 | // starguider data points:
|
---|
33 | //
|
---|
34 | // * MPointingDevCalc.NumMinStars: 8
|
---|
35 | // Minimum number of identified stars required to accep the data
|
---|
36 | //
|
---|
37 | // * MPointingDevCalc.NsbLevel: 3.0
|
---|
38 | // Minimum deviation (in rms) from the the mean allowed for the
|
---|
39 | // measured NSB (noise in the ccd camera)
|
---|
40 | //
|
---|
41 | // * MPointingDevCalc.NsbMin: 30
|
---|
42 | // - minimum NSB to be used in mean/rms calculation
|
---|
43 | //
|
---|
44 | // * MPointingDevCalc.NsbMax: 60
|
---|
45 | // - maximum NSB to be used in mean/rms calculation
|
---|
46 | //
|
---|
47 | // * MPointingDevCalc.MaxAbsDev: 15
|
---|
48 | // - Maximum absolute deviation which is consideres as valid (arcmin)
|
---|
49 | //
|
---|
50 | // Starguider data which doens't fullfill this requirements is ignored.
|
---|
51 | // If the measures NSB==0 (file too old, starguider didn't write down
|
---|
52 | // these values) the checks based on NSB and NumMinStar are skipped.
|
---|
53 | //
|
---|
54 | // The calculation of NSB mean and rms is reset for each file (ReInit)
|
---|
55 | //
|
---|
56 | // If your starguider data doesn't fullfill this requirement the latest
|
---|
57 | // value which could be correctly calculated is used instead. If the time
|
---|
58 | // for which no valid value can be calculated exceeds one minute
|
---|
59 | // the return value is reset to 0/0. The maximum time allowed without
|
---|
60 | // a valid value can be setup using:
|
---|
61 | //
|
---|
62 | // * MPointingDevCalc.MaxAge: 1
|
---|
63 | // Maximum time before the starguider is reset to 0/0 in minutes
|
---|
64 | //
|
---|
65 | // Note, that the starguider itself is not well calibrated. Therefore
|
---|
66 | // it is necessary to do a starguider calibration in our software.
|
---|
67 | //
|
---|
68 | // There are two options:
|
---|
69 | //
|
---|
70 | // * Simple starguider calibration using offsets in the camera plane
|
---|
71 | //
|
---|
72 | // The starguider is calibrated by taking its values (dZd/dAz)
|
---|
73 | // adding them to the source position, calculating the source position
|
---|
74 | // in the camera plane and adding the offsets. To switch off the
|
---|
75 | // full starguider calibration do:
|
---|
76 | //
|
---|
77 | // * MPointingDevCalc.PointingModels:
|
---|
78 | //
|
---|
79 | // To set the offsets (in units of degree) use
|
---|
80 | //
|
---|
81 | // * MPointingDevCalc.Dx: -0.001
|
---|
82 | // * MPointingDevCalc.Dy: -0.004
|
---|
83 | //
|
---|
84 | // * A starguider calibration using a pointing model calculated
|
---|
85 | // from calibration data, so called TPoints
|
---|
86 | //
|
---|
87 | // Because the pointing model can change from time to time
|
---|
88 | // you can give the run-number from which on a new pointing
|
---|
89 | // model is valid. The run itself is included, e.g.:
|
---|
90 | //
|
---|
91 | // * MPointingDevCalc.PointingModels: 85240 89180
|
---|
92 | // * MPointingDevCalc.FilePrefix: resources/starguider
|
---|
93 | //
|
---|
94 | // mean that for all runs<85240 the simple offset correction is used.
|
---|
95 | // For runs >=85240 and <89180 the file resources/starguider0085240.txt
|
---|
96 | // and for runs >=89180 the file resources/starguider0089180.txt is
|
---|
97 | // used. To setup a default file for all runs before 85240 setup
|
---|
98 | // a low number (eg. 0 or 1)
|
---|
99 | //
|
---|
100 | // In the case a pointing model is used additional offsets in
|
---|
101 | // the x/y-camera plane (in units of deg) can be set using the DX
|
---|
102 | // and DY parameters of the pointing model. The fDx and fDy data
|
---|
103 | // members of this class are ignored. To overwrite the starguider
|
---|
104 | // calibrated offset in either Az or Zd with a constant, you
|
---|
105 | // can use the PX/PY directive in the pointing model. (To enable
|
---|
106 | // the overwrite set the third column, the error, to a value
|
---|
107 | // greater than zero)
|
---|
108 | //
|
---|
109 | //
|
---|
110 | // At the PostProcessing step a table with statistics is print if the
|
---|
111 | // debug level is greater or equal 3 (in most applications it is switched
|
---|
112 | // on by -v3)
|
---|
113 | //
|
---|
114 | //
|
---|
115 | // Pointing Models:
|
---|
116 | // ----------------
|
---|
117 | //
|
---|
118 | // What we know so far about (maybe) important changes in cosy or to the
|
---|
119 | // hardware:
|
---|
120 | //
|
---|
121 | // 25.01.2010: Changed scaling factors for TPoint camera (M1)
|
---|
122 | //
|
---|
123 | // 11.01.2010: New TPoint camera (M1)
|
---|
124 | //
|
---|
125 | // 18.03.2006: The camera holding has been repaired and the camera got
|
---|
126 | // shifted a little bit.
|
---|
127 | //
|
---|
128 | // 16.04.2006: Around this date a roque lamp focussing has been done
|
---|
129 | //
|
---|
130 | // 25.04.2006: A missalignment introduced with the roque adjust has been
|
---|
131 | // corrected
|
---|
132 | //
|
---|
133 | // The starguider pointing model for the time before 18.3.2006 and after
|
---|
134 | // April 2006 (in fact there are no TPoints until 07/2006 to check it)
|
---|
135 | // and for the period 07/2006 to (at least) 06/2007 are very similar.
|
---|
136 | //
|
---|
137 | // The pointing model for the time between 18.3.2006 and 04/2006 is
|
---|
138 | // clearly different, mainly giving different Azimuth values between
|
---|
139 | // Zenith and roughly ~25deg, and a slight offset on both axes.
|
---|
140 | //
|
---|
141 | // 10.5.2006: pos1 -= pos0 commented (What was the mentioned fix?)
|
---|
142 | // 29.6.2006: repaired
|
---|
143 | //
|
---|
144 | // 23.3.2006: new starguider algorithm
|
---|
145 | //
|
---|
146 | // 17.3.2005: Fixed units of "nompos" in MDriveCom
|
---|
147 | //
|
---|
148 | //
|
---|
149 | // New pointing models have been installed (if the pointing model
|
---|
150 | // is different, than the previous one it might mean, that also
|
---|
151 | // the starguider calibration is different.) The date only means
|
---|
152 | // day-time (noon) at which the model has been changed.
|
---|
153 | //
|
---|
154 | // 29. Apr. 2004 ~25800
|
---|
155 | // 5. Aug. 2004 ~32000
|
---|
156 | // 19. Aug. 2004 ~33530
|
---|
157 | // 7. Jun. 2005 ~57650
|
---|
158 | // 8. Jun. 2005
|
---|
159 | // 9. Jun. 2005 ~57860
|
---|
160 | // 12. Sep. 2005 ~68338
|
---|
161 | // 24. Nov. 2005 ~75562
|
---|
162 | // 17. Oct. 2006 ~103130
|
---|
163 | // 17. Jun. 2007 ~248193
|
---|
164 | // 18. Oct. 2007 291104 // Correction for the offsets introduced by AMC
|
---|
165 | // 14. Jan. 2008 328198 // Complete new pointing model
|
---|
166 | // 11. Jun. 2008 1000534 (ca. 23:00) // Before new TPoints
|
---|
167 | // 19. Jun. 2008 (ca. 15:00) // From new TPoints
|
---|
168 | // 7. Mar. 2009 (ca. 14:00) // From new TPoints (0808-0902)
|
---|
169 | // 14. May 2009 // M1/M2 after upgrade (from TPoints taken in the days before)
|
---|
170 | // 17. Aug. 2009 New pointing models for both telescopes
|
---|
171 | // 01. Feb. 2010 New pointing models for both telescopes
|
---|
172 | // 26. Feb. 2010 New pointing models for both telescopes
|
---|
173 | // 31. Mar. 2010 New pointing models for both telescopes
|
---|
174 | // 29. Sep. 2010 New pointing model MAGIC I
|
---|
175 | // 01. Dec. 2010 New pointing models for both telescopes
|
---|
176 | // 06. Jan. 2011 New pointing models for both telescopes
|
---|
177 | //
|
---|
178 | // From 2.2.2006 beginnig of the night (21:05h, run >=81855) to 24.2.2006
|
---|
179 | // beginning of the the night (20:34h, run<83722) the LEDs did not work.
|
---|
180 | // In the time after this incident the shift crew often forgot to switch on
|
---|
181 | // the LEDs at the beginning of the night!
|
---|
182 | //
|
---|
183 | // [2006-03-07 00:10:58] In the daytime, we raised the position of the
|
---|
184 | // 9 o'clock LED by one screw hole to make it visible when the TPoint
|
---|
185 | // Lid is closed. (< run 84980)
|
---|
186 | //
|
---|
187 | // Mirror refocussing around 23.Apr.2006, from the Runbook.
|
---|
188 | //
|
---|
189 | // 25./26.4.2006: Run 89180
|
---|
190 | //
|
---|
191 | // (Note: The year here is not a typo!)
|
---|
192 | // [2007-04-25 23:50:39]
|
---|
193 | // Markus is performing AMC focussing.
|
---|
194 | //
|
---|
195 | // Mirror refocussing around 4.Aug.2007, from the Runbook:
|
---|
196 | //
|
---|
197 | // [2007-08-04 04:46:47]
|
---|
198 | // We finished with the focussing with Polaris. The images need to be
|
---|
199 | // analysed and new LUTs generated.
|
---|
200 | //
|
---|
201 | // [2007-08-04 23:47:30]
|
---|
202 | // Actually we see that the mispointing is always large; probably since
|
---|
203 | // the LUT tables have not yet been adjusted to the new focussing.
|
---|
204 | //
|
---|
205 | // [2007-08-03 23:07:58]
|
---|
206 | // Data taking stopped. Mirror focussing.
|
---|
207 | //
|
---|
208 | // [2007-08-05 00:09:16]
|
---|
209 | // We take some pictures on stars nearby Cyg X3 with the sbig camera;
|
---|
210 | // actually the spot doesn't look very nice... The pictures have been
|
---|
211 | // saved with name Deneb- and Sadr- Polaris seems a bit better. Should we
|
---|
212 | // have new LUT tables after the focussing?
|
---|
213 | //
|
---|
214 | // [2007-08-10 20:18:48]
|
---|
215 | // Tonight we take first images of Polaris with a new LUT file generated
|
---|
216 | // based on the recent focussing. The image will be analysed tomorrow and
|
---|
217 | // than new LUTs will be generated. For tonight the focussing is still not
|
---|
218 | // changed.
|
---|
219 | //
|
---|
220 | // [2007-08-14 20:57:59]
|
---|
221 | // The weather is fine. There is a group of hobby astronomers at the
|
---|
222 | // helicopter parking. Before data taking, we tried to check the new LUTs.
|
---|
223 | // However, because of technical problems with the new LUTs we had to
|
---|
224 | // postpone the measurements. We lost some 10 min of data taking because
|
---|
225 | // of this.
|
---|
226 | //
|
---|
227 | // [2007-08-14 22:29:37] Run 267253
|
---|
228 | // Before continuing the observation we perform a focussing test with the
|
---|
229 | // new LUTs from recent Polaris focussing. Note: Data on Her X-1 was taken
|
---|
230 | // with old focussing. We performed PSF measurements on Kornephorus (Zd
|
---|
231 | // 31.77, Az 264,67) first with old LUTs and then with new LUTs. We took
|
---|
232 | // T-points with both focussing. The first T-Point corresponds to the
|
---|
233 | // previous focussing and the second with the improved. Please check both
|
---|
234 | // T-points for eventual misspointing. We found big improvement of the PSF
|
---|
235 | // with the new LUTs and will therefore continue from now on with the new
|
---|
236 | // focussing. Having a first look at the SBIG pictures we see a slightly
|
---|
237 | // misspointing of ~0.1 deg with the new LUTs.
|
---|
238 | //
|
---|
239 | // [2007-08-14 22:46:10]
|
---|
240 | // Comparing the trigger rate with yesterday night we do not see an
|
---|
241 | // improvement with the new focussing.
|
---|
242 | //
|
---|
243 | // [2007-10-27 email]
|
---|
244 | // [...]
|
---|
245 | // When removing the cover of the motor all 4 nuts which fix the gear of
|
---|
246 | // the motor to the structure fell out. The motor was completely loose.
|
---|
247 | // [...]
|
---|
248 | // We checked also the second azimuth motor and .. the same situation.
|
---|
249 | // [...]
|
---|
250 | // Peter Sawallisch opened and fixed all 3 motors of MAGIC-I. We checked
|
---|
251 | // all critical screws that came to our minds and fixed them as much as
|
---|
252 | // possible with loctite and counter nuts. No damage whatsoever has been
|
---|
253 | // found.
|
---|
254 | //
|
---|
255 | // [2009-08-05 Mail from Markus Garcz. about M1]
|
---|
256 | // [...]
|
---|
257 | // On the 23.07.2009 we tested the new LUTs using the star Etamin. The
|
---|
258 | // PSF (sigma) improved from 11.0 mm to 9.5 mm and the total light content
|
---|
259 | // in one PMT from 54% to 62%.
|
---|
260 | // Furthermore I saw a small movement of the spot after the new focussing:
|
---|
261 | // dX = 3 CCD pixel = 7 mm
|
---|
262 | // dY = 2 CCD pixel = 4.7 mm
|
---|
263 | // The new LUTs are now installed since the 23.07.2009 and are used as
|
---|
264 | // default during the observation.
|
---|
265 | // [...]
|
---|
266 | //
|
---|
267 | //
|
---|
268 | // [2009-10-16 Email Adrian]
|
---|
269 | //
|
---|
270 | // I checked the AMC2 log files to check which LUTs had been used in which nights
|
---|
271 | // until 09/09/12 the old version from March had always been used
|
---|
272 | // 09/09/13-09/09/16 LUT_090913 (test) was used (by mistake)
|
---|
273 | // 09/09/17-09/09/20 switched back to the version from March
|
---|
274 | // since 09/09/21 using actual version LUT_090918
|
---|
275 | // (PSF at ~10deg still rather poor, but I see strange
|
---|
276 | // effects I do not yet understand and therefore cannot
|
---|
277 | // correct)
|
---|
278 | // Especially: nothing has changed in the AMC on 09/09/09,
|
---|
279 | // except that around that date the SBIG camera was readjusted (but this
|
---|
280 | // has no direct affect on the AMC)
|
---|
281 | //
|
---|
282 | //
|
---|
283 | // [2010-02-03 M. Garcz.]
|
---|
284 | //
|
---|
285 | // New LUTs for M2
|
---|
286 | //
|
---|
287 | // [2010-06-14 M. Garcz.]
|
---|
288 | //
|
---|
289 | // New LUTs for M1
|
---|
290 | //
|
---|
291 | // [2010-07-02 M. Garcz.]
|
---|
292 | //
|
---|
293 | // New LUTs for M1
|
---|
294 | //
|
---|
295 | // [2010-10-24]
|
---|
296 | //
|
---|
297 | // Since about this date we see a higher spread in the TPoints. The reason we
|
---|
298 | // found is that SA sends the coordinates of the orginal position and not the
|
---|
299 | // target position to the AMC thus the mirrors are not well focussed.
|
---|
300 | // This should mainly effect the quality of the TPoints but not the pointing
|
---|
301 | // model. However, M2 shows a slight offset.
|
---|
302 | //
|
---|
303 | // [2010-11-11 M. Garcz.]
|
---|
304 | //
|
---|
305 | // Last night, TPoints for M2 taken without target
|
---|
306 | //
|
---|
307 | // [2010-11-18 M. Garcz.]
|
---|
308 | //
|
---|
309 | // New LUTs for M1
|
---|
310 | //
|
---|
311 | //
|
---|
312 | // Others
|
---|
313 | // ------
|
---|
314 | //
|
---|
315 | // The pointing of both(!) telescopes changed >2009/06/11 2:36h
|
---|
316 | // In both cases the reason is unknown.
|
---|
317 | //
|
---|
318 | // The pointing has changed for M2 after 2009/09/09 noon
|
---|
319 | // It recovered to the previous pointing at 2009/10/08 noon
|
---|
320 | // Resons unknown (there were LUT changes but at different dates)
|
---|
321 | //
|
---|
322 | //
|
---|
323 | // ToDo:
|
---|
324 | // -----
|
---|
325 | //
|
---|
326 | // * Is 0/0 the best assumption if the starguider partly fails?
|
---|
327 | //
|
---|
328 | //
|
---|
329 | // Input Container:
|
---|
330 | // MRawRunHeader
|
---|
331 | // MReportStarguider
|
---|
332 | //
|
---|
333 | // Output Container:
|
---|
334 | // MPointingDev
|
---|
335 | //
|
---|
336 | /////////////////////////////////////////////////////////////////////////////
|
---|
337 | #include "MPointingDevCalc.h"
|
---|
338 |
|
---|
339 | #include <stdlib.h> // atoi (Ubuntu 8.10)
|
---|
340 |
|
---|
341 | #include "MLog.h"
|
---|
342 | #include "MLogManip.h"
|
---|
343 |
|
---|
344 | #include "MParList.h"
|
---|
345 |
|
---|
346 | #include "MAstro.h"
|
---|
347 | #include "MPointing.h"
|
---|
348 | #include "MPointingDev.h"
|
---|
349 | #include "MRawRunHeader.h"
|
---|
350 | #include "MReportStarguider.h"
|
---|
351 |
|
---|
352 | ClassImp(MPointingDevCalc);
|
---|
353 |
|
---|
354 | using namespace std;
|
---|
355 |
|
---|
356 | const TString MPointingDevCalc::fgFilePrefix="resources/starguider";
|
---|
357 |
|
---|
358 | // --------------------------------------------------------------------------
|
---|
359 | //
|
---|
360 | // Destructor. Call Clear() and delete fPointingModels if any.
|
---|
361 | //
|
---|
362 | MPointingDevCalc::~MPointingDevCalc()
|
---|
363 | {
|
---|
364 | Clear();
|
---|
365 | }
|
---|
366 |
|
---|
367 | // --------------------------------------------------------------------------
|
---|
368 | //
|
---|
369 | // Delete fPointing and set pointing to NULL
|
---|
370 | //
|
---|
371 | void MPointingDevCalc::Clear(Option_t *o)
|
---|
372 | {
|
---|
373 | if (fPointing)
|
---|
374 | delete fPointing;
|
---|
375 |
|
---|
376 | fPointing = NULL;
|
---|
377 | }
|
---|
378 |
|
---|
379 | // --------------------------------------------------------------------------
|
---|
380 | //
|
---|
381 | // Sort the entries in fPoinitngModels
|
---|
382 | //
|
---|
383 | void MPointingDevCalc::SortPointingModels()
|
---|
384 | {
|
---|
385 | const int n = fPointingModels.GetSize();
|
---|
386 |
|
---|
387 | TArrayI idx(n);
|
---|
388 |
|
---|
389 | TMath::Sort(n, fPointingModels.GetArray(), idx.GetArray(), kFALSE);
|
---|
390 |
|
---|
391 | const TArrayI arr(fPointingModels);
|
---|
392 |
|
---|
393 | for (int i=0; i<n; i++)
|
---|
394 | fPointingModels[i] = arr[idx[i]];
|
---|
395 | }
|
---|
396 |
|
---|
397 | // --------------------------------------------------------------------------
|
---|
398 | //
|
---|
399 | // Set new pointing models
|
---|
400 | //
|
---|
401 | void MPointingDevCalc::SetPointingModels(const TString &models)
|
---|
402 | {
|
---|
403 | fPointingModels.Set(0);
|
---|
404 |
|
---|
405 | if (models.IsNull())
|
---|
406 | return;
|
---|
407 |
|
---|
408 | TObjArray *arr = models.Tokenize(" ");
|
---|
409 |
|
---|
410 | const int n = arr->GetEntries();
|
---|
411 | fPointingModels.Set(n);
|
---|
412 |
|
---|
413 | for (int i=0; i<n; i++)
|
---|
414 | fPointingModels[i] = atoi((*arr)[i]->GetName());
|
---|
415 |
|
---|
416 | delete arr;
|
---|
417 |
|
---|
418 | SortPointingModels();
|
---|
419 | }
|
---|
420 |
|
---|
421 | // --------------------------------------------------------------------------
|
---|
422 | //
|
---|
423 | // Return a string with the pointing models, seperated by a space.
|
---|
424 | //
|
---|
425 | TString MPointingDevCalc::GetPointingModels() const
|
---|
426 | {
|
---|
427 | TString rc;
|
---|
428 | for (int i=0; i<fPointingModels.GetSize(); i++)
|
---|
429 | rc += Form ("%d ", fPointingModels[i]);
|
---|
430 |
|
---|
431 | return rc;
|
---|
432 | }
|
---|
433 |
|
---|
434 | // --------------------------------------------------------------------------
|
---|
435 | //
|
---|
436 | // Add a number to the pointing models
|
---|
437 | //
|
---|
438 | void MPointingDevCalc::AddPointingModel(UInt_t runnum)
|
---|
439 | {
|
---|
440 | const int n = fPointingModels.GetSize();
|
---|
441 | for (int i=0; i<n; i++)
|
---|
442 | if ((UInt_t)fPointingModels[i]==runnum)
|
---|
443 | {
|
---|
444 | *fLog << warn << "WARNING - Pointing model " << runnum << " already in list... ignored." << endl;
|
---|
445 | return;
|
---|
446 | }
|
---|
447 |
|
---|
448 | fPointingModels.Set(n+1);
|
---|
449 | fPointingModels[n] = runnum;
|
---|
450 |
|
---|
451 | SortPointingModels();
|
---|
452 | }
|
---|
453 |
|
---|
454 | // --------------------------------------------------------------------------
|
---|
455 | //
|
---|
456 | // Find the highest number in the array which is lower or equal num.
|
---|
457 | //
|
---|
458 | UInt_t MPointingDevCalc::FindPointingModel(UInt_t num)
|
---|
459 | {
|
---|
460 | const int n = fPointingModels.GetSize();
|
---|
461 | if (n==0)
|
---|
462 | return (UInt_t)-1;
|
---|
463 |
|
---|
464 | // Loop over all pointing models
|
---|
465 | for (int i=0; i<n; i++)
|
---|
466 | {
|
---|
467 | // The number stored in the array did not yet overtake the runnumber
|
---|
468 | if ((UInt_t)fPointingModels[i]<num)
|
---|
469 | continue;
|
---|
470 |
|
---|
471 | // The first pointing model is later than this run: use a default
|
---|
472 | if (i==0)
|
---|
473 | return (UInt_t)-1;
|
---|
474 |
|
---|
475 | // The last entry in the array is the right one.
|
---|
476 | return fPointingModels[i-1];
|
---|
477 | }
|
---|
478 |
|
---|
479 | // Runnumber is after last entry of pointing models. Use the last one.
|
---|
480 | return fPointingModels[n-1];
|
---|
481 | }
|
---|
482 |
|
---|
483 | // --------------------------------------------------------------------------
|
---|
484 | //
|
---|
485 | // Clear the pointing model. If run-number >= 87751 read the new
|
---|
486 | // pointing model with fFilePrefix
|
---|
487 | //
|
---|
488 | Bool_t MPointingDevCalc::ReadPointingModel(const MRawRunHeader &run)
|
---|
489 | {
|
---|
490 | const UInt_t num = FindPointingModel(run.GetRunNumber());
|
---|
491 |
|
---|
492 | // No poinitng models are defined. Use simple dx/dy-calibration
|
---|
493 | if (num==(UInt_t)-1)
|
---|
494 | {
|
---|
495 | Clear();
|
---|
496 | return kTRUE;
|
---|
497 | }
|
---|
498 |
|
---|
499 | // compile the name for the starguider files
|
---|
500 | // The file with the number 00000000 is the default file
|
---|
501 | TString fname = Form("%s%08d.txt", fFilePrefix.Data(), num);
|
---|
502 |
|
---|
503 | if (!fPointing)
|
---|
504 | fPointing = new MPointing;
|
---|
505 |
|
---|
506 | if (fname==fPointing->GetName())
|
---|
507 | {
|
---|
508 | *fLog << inf << fname << " already loaded." << endl;
|
---|
509 | return kTRUE;
|
---|
510 | }
|
---|
511 |
|
---|
512 | return fPointing->Load(fname);
|
---|
513 | }
|
---|
514 |
|
---|
515 | // --------------------------------------------------------------------------
|
---|
516 | //
|
---|
517 | // Check the file/run type from the run-header and if it is a data file
|
---|
518 | // load starguider calibration.
|
---|
519 | //
|
---|
520 | Bool_t MPointingDevCalc::ReInit(MParList *plist)
|
---|
521 | {
|
---|
522 | MRawRunHeader *run = (MRawRunHeader*)plist->FindObject("MRawRunHeader");
|
---|
523 | if (!run)
|
---|
524 | {
|
---|
525 | *fLog << err << "MRawRunHeader not found... aborting." << endl;
|
---|
526 | return kFALSE;
|
---|
527 | }
|
---|
528 |
|
---|
529 | fNsbSum = 0;
|
---|
530 | fNsbSq = 0;
|
---|
531 | fNsbCount = 0;
|
---|
532 |
|
---|
533 | fRunType = run->GetRunType();
|
---|
534 |
|
---|
535 | switch (fRunType)
|
---|
536 | {
|
---|
537 | case MRawRunHeader::kRTData:
|
---|
538 | if (!fReport)
|
---|
539 | *fLog << warn << "MReportStarguider not found... skipped." << endl;
|
---|
540 | return ReadPointingModel(*run);
|
---|
541 |
|
---|
542 | case MRawRunHeader::kRTMonteCarlo:
|
---|
543 | return kTRUE;
|
---|
544 |
|
---|
545 | case MRawRunHeader::kRTPedestal:
|
---|
546 | *fLog << err << "Cannot work in a pedestal Run!... aborting." << endl;
|
---|
547 | return kFALSE;
|
---|
548 |
|
---|
549 | case MRawRunHeader::kRTCalibration:
|
---|
550 | *fLog << err << "Cannot work in a calibration Run!... aborting." << endl;
|
---|
551 | return kFALSE;
|
---|
552 |
|
---|
553 | default:
|
---|
554 | *fLog << err << "Run Type " << fRunType << " unknown!... aborting." << endl;
|
---|
555 | return kFALSE;
|
---|
556 | }
|
---|
557 |
|
---|
558 | return kTRUE;
|
---|
559 | }
|
---|
560 |
|
---|
561 | // --------------------------------------------------------------------------
|
---|
562 | //
|
---|
563 | // Search for 'MPointingPos'. Create if not found.
|
---|
564 | //
|
---|
565 | Int_t MPointingDevCalc::PreProcess(MParList *plist)
|
---|
566 | {
|
---|
567 | fDeviation = (MPointingDev*)plist->FindCreateObj("MPointingDev");
|
---|
568 | fReport = (MReportStarguider*)plist->FindObject("MReportStarguider");
|
---|
569 |
|
---|
570 | // We use kRTNone here as a placeholder for data runs.
|
---|
571 | fRunType = MRawRunHeader::kRTNone;
|
---|
572 | fLastMjd = -1;
|
---|
573 |
|
---|
574 | fSkip.Reset();
|
---|
575 |
|
---|
576 | // In cases in which ReInit isn't called right in time (e.g. if
|
---|
577 | // the first starguider event comes before the first data event)
|
---|
578 | fNsbSum = 0;
|
---|
579 | fNsbSq = 0;
|
---|
580 | fNsbCount = 0;
|
---|
581 |
|
---|
582 | return fDeviation ? kTRUE : kFALSE;
|
---|
583 | }
|
---|
584 |
|
---|
585 | // --------------------------------------------------------------------------
|
---|
586 | //
|
---|
587 | // Increase fSkip[i] by one. If the data in fDeviation is outdated (older
|
---|
588 | // than fMaxAge) and the current report should be skipped reset DevZdAz and
|
---|
589 | // DevXY and fSkip[6] is increased by one.
|
---|
590 | //
|
---|
591 | void MPointingDevCalc::Skip(Int_t i)
|
---|
592 | {
|
---|
593 | fSkip[i]++;
|
---|
594 |
|
---|
595 | const Double_t diff = (fReport->GetMjd()-fLastMjd)*1440; // [min] 1440=24*60
|
---|
596 | if (diff<fMaxAge && fLastMjd>0)
|
---|
597 | return;
|
---|
598 |
|
---|
599 | fDeviation->SetDevZdAz(0, 0);
|
---|
600 | fDeviation->SetDevXY(0, 0);
|
---|
601 | fSkip[6]++;
|
---|
602 | }
|
---|
603 |
|
---|
604 | // --------------------------------------------------------------------------
|
---|
605 | //
|
---|
606 | // Do a full starguider calibration using a pointing model for the starguider.
|
---|
607 | //
|
---|
608 | void MPointingDevCalc::DoCalibration(Double_t devzd, Double_t devaz) const
|
---|
609 | {
|
---|
610 | if (!fPointing)
|
---|
611 | {
|
---|
612 | // Do a simple starguider calibration using a simple offset in x and y
|
---|
613 | fDeviation->SetDevZdAz(devzd, devaz);
|
---|
614 |
|
---|
615 | // Linear starguider calibration taken from April/May data
|
---|
616 | // For calibration add MDriveReport::GetErrorZd/Az !
|
---|
617 | fDeviation->SetDevXY(fDx, fDy); // 1arcmin ~ 5mm
|
---|
618 |
|
---|
619 | return;
|
---|
620 | }
|
---|
621 |
|
---|
622 | // Get the nominal position the star is at the sky
|
---|
623 | // Unit: deg
|
---|
624 | ZdAz nom(fReport->GetNominalZd(), fReport->GetNominalAz());
|
---|
625 | nom *= TMath::DegToRad();
|
---|
626 |
|
---|
627 | // Get the mispointing measured by the telescope. It is
|
---|
628 | // calculate as follows:
|
---|
629 | //
|
---|
630 | // The mispointing measured by the starguider:
|
---|
631 | // ZdAz mis(devzd, devaz);
|
---|
632 | // mis *= TMath::DegToRad();
|
---|
633 |
|
---|
634 | // The pointing model is the conversion from the real pointing
|
---|
635 | // position of the telescope into the pointing position measured
|
---|
636 | // by the starguider.
|
---|
637 | //
|
---|
638 | // To keep as close to the fitted model we use the forward correction.
|
---|
639 |
|
---|
640 | // Position at which the starguider camera is pointing in real:
|
---|
641 | // pointing position = nominal position - dev
|
---|
642 | //
|
---|
643 | // The position measured as the starguider's pointing position
|
---|
644 | ZdAz pos(nom); // cpos = sao - dev
|
---|
645 | pos -= ZdAz(devzd, devaz)*TMath::DegToRad();
|
---|
646 |
|
---|
647 | // Now we convert the starguider's pointing position into the
|
---|
648 | // telescope pointing position (the pointing model converts
|
---|
649 | // the telescope pointing position into the starguider pointing
|
---|
650 | // position)
|
---|
651 | ZdAz point = fPointing->CorrectBack(pos); //FWD!!!
|
---|
652 |
|
---|
653 | // MSrcPosCalc uses the following condition to calculate the
|
---|
654 | // source position in the camera:
|
---|
655 | // real pointing pos = nominal pointing pos - dev
|
---|
656 | // --> dev = nominal - real
|
---|
657 | // Therefor we calculate dev as follows:
|
---|
658 | ZdAz dev(nom);
|
---|
659 | dev -= point;
|
---|
660 | dev *= TMath::RadToDeg();
|
---|
661 |
|
---|
662 | /*
|
---|
663 | // We chose the other way. It is less accurate because is is the
|
---|
664 | // other was than the poinitng model was fittet, but it is more
|
---|
665 | // accurate because the nominal (i.e. real) pointing position
|
---|
666 | // is less accurately known than the position returned by the
|
---|
667 | // starguider.
|
---|
668 | //
|
---|
669 | // Calculate the deviation which would be measured by the starguider
|
---|
670 | // if applied to a perfectly pointing telescope.
|
---|
671 | ZdAz dev = fPointing->Correct(nom);
|
---|
672 | dev -= nom;
|
---|
673 |
|
---|
674 | // Now add these offsets and the starguider measured offsets to
|
---|
675 | // the real pointing deviation of the telescope (note, that
|
---|
676 | // signs here are just conventions)
|
---|
677 | dev += ZdAz(devzd, devaz)*TMath::DegToRad(); // --> nom-mis
|
---|
678 | dev *= TMath::RadToDeg();
|
---|
679 | */
|
---|
680 |
|
---|
681 | // Check if the starguider pointing model requests overwriting
|
---|
682 | // of the values with constants (e.g. 0)
|
---|
683 | devaz = fPointing->IsPxValid() ? fPointing->GetPx() : dev.Az();
|
---|
684 | devzd = fPointing->IsPyValid() ? fPointing->GetPy() : dev.Zd();
|
---|
685 |
|
---|
686 | fDeviation->SetDevZdAz(devzd, devaz);
|
---|
687 | fDeviation->SetDevXY(fPointing->GetDxy());
|
---|
688 | }
|
---|
689 |
|
---|
690 | Int_t MPointingDevCalc::ProcessStarguiderReport()
|
---|
691 | {
|
---|
692 | Double_t devzd = fReport->GetDevZd(); // [arcmin]
|
---|
693 | Double_t devaz = fReport->GetDevAz(); // [arcmin]
|
---|
694 | if (devzd==0 && devaz==0)
|
---|
695 | {
|
---|
696 | Skip(1);
|
---|
697 | return kTRUE;
|
---|
698 | }
|
---|
699 |
|
---|
700 | if (!fReport->IsMonitoring())
|
---|
701 | {
|
---|
702 | Skip(2);
|
---|
703 | return kTRUE;
|
---|
704 | }
|
---|
705 |
|
---|
706 | devzd /= 60; // Convert from arcmin to deg
|
---|
707 | devaz /= 60; // Convert from arcmin to deg
|
---|
708 |
|
---|
709 | const Double_t nsb = fReport->GetSkyBrightness();
|
---|
710 | if (nsb>0)
|
---|
711 | {
|
---|
712 | if (nsb>fNsbMin && nsb<fNsbMax)
|
---|
713 | {
|
---|
714 | fNsbSum += nsb;
|
---|
715 | fNsbSq += nsb*nsb;
|
---|
716 | fNsbCount++;
|
---|
717 | }
|
---|
718 |
|
---|
719 | if (fNsbCount>0)
|
---|
720 | {
|
---|
721 | const Double_t sum = fNsbSum/fNsbCount;
|
---|
722 | const Double_t sq = fNsbSq /fNsbCount;
|
---|
723 |
|
---|
724 | const Double_t rms = fNsbLevel*TMath::Sqrt(sq - sum*sum);
|
---|
725 |
|
---|
726 | if (nsb<sum-rms || nsb>sum+rms)
|
---|
727 | {
|
---|
728 | Skip(3);
|
---|
729 | return kTRUE;
|
---|
730 | }
|
---|
731 | }
|
---|
732 |
|
---|
733 | if (fReport->GetNumIdentifiedStars()<fNumMinStars)
|
---|
734 | {
|
---|
735 | Skip(4);
|
---|
736 | return kTRUE;
|
---|
737 | }
|
---|
738 | }
|
---|
739 |
|
---|
740 | // >= 87751 (31.3.06 12:00)
|
---|
741 |
|
---|
742 | // Calculate absolute deviation
|
---|
743 | const Double_t dev = MAstro::GetDevAbs(fReport->GetNominalZd(),
|
---|
744 | fReport->GetNominalZd()-fReport->GetDevZd()/60,
|
---|
745 | devaz);
|
---|
746 |
|
---|
747 | // Sanity check... larger deviation are strange and ignored
|
---|
748 | if (dev*60>fMaxAbsDev)
|
---|
749 | {
|
---|
750 | Skip(5);
|
---|
751 | return kTRUE;
|
---|
752 | }
|
---|
753 |
|
---|
754 | DoCalibration(devzd, devaz);
|
---|
755 |
|
---|
756 | fSkip[0]++;
|
---|
757 | fLastMjd = fReport->GetMjd();
|
---|
758 |
|
---|
759 | return kTRUE;
|
---|
760 | }
|
---|
761 |
|
---|
762 | // --------------------------------------------------------------------------
|
---|
763 | //
|
---|
764 | // See class description.
|
---|
765 | //
|
---|
766 | Int_t MPointingDevCalc::Process()
|
---|
767 | {
|
---|
768 | switch (fRunType)
|
---|
769 | {
|
---|
770 | case MRawRunHeader::kRTNone:
|
---|
771 | case MRawRunHeader::kRTData:
|
---|
772 | return fReport ? ProcessStarguiderReport() : kTRUE;
|
---|
773 |
|
---|
774 | case MRawRunHeader::kRTMonteCarlo:
|
---|
775 | fSkip[0]++;
|
---|
776 | fDeviation->SetDevZdAz(0, 0);
|
---|
777 | fDeviation->SetDevXY(0, 0);
|
---|
778 | return kTRUE;
|
---|
779 | }
|
---|
780 | return kTRUE;
|
---|
781 | }
|
---|
782 |
|
---|
783 | // --------------------------------------------------------------------------
|
---|
784 | //
|
---|
785 | // Print execution statistics
|
---|
786 | //
|
---|
787 | Int_t MPointingDevCalc::PostProcess()
|
---|
788 | {
|
---|
789 | if (GetNumExecutions()==0)
|
---|
790 | return kTRUE;
|
---|
791 |
|
---|
792 | *fLog << inf << endl;
|
---|
793 | *fLog << GetDescriptor() << " execution statistics:" << endl;
|
---|
794 | PrintSkipped(fSkip[1], "Starguider deviation not set, is exactly 0/0");
|
---|
795 | PrintSkipped(fSkip[2], "Starguider was not monitoring (eg. LEDs off)");
|
---|
796 | PrintSkipped(fSkip[3], Form("NSB out of %.1f sigma range", fNsbLevel));
|
---|
797 | PrintSkipped(fSkip[4], Form("Number of identified stars < %d", fNumMinStars));
|
---|
798 | PrintSkipped(fSkip[5], Form("Absolute deviation > %.1farcmin", fMaxAbsDev));
|
---|
799 | PrintSkipped(fSkip[6], Form("Events set to 0 because older than %.1fmin", fMaxAge));
|
---|
800 | *fLog << " " << (int)fSkip[0] << " (" << Form("%5.1f", 100.*fSkip[0]/GetNumExecutions()) << "%) Evts survived calculation!" << endl;
|
---|
801 | *fLog << endl;
|
---|
802 |
|
---|
803 | return kTRUE;
|
---|
804 | }
|
---|
805 |
|
---|
806 | // --------------------------------------------------------------------------
|
---|
807 | //
|
---|
808 | // MPointingDevCalc.NumMinStars: 8
|
---|
809 | // MPointingDevCalc.NsbLevel: 3.0
|
---|
810 | // MPointingDevCalc.NsbMin: 30
|
---|
811 | // MPointingDevCalc.NsbMax: 60
|
---|
812 | // MPointingDevCalc.MaxAbsDev: 15
|
---|
813 | // MPointingDevCalc.MaxAge: 1.0
|
---|
814 | // MPointingDevCalc.Dx: -0.001
|
---|
815 | // MPointingDevCalc.Dy: -0.004
|
---|
816 | //
|
---|
817 | // For a detailed description see the class reference.
|
---|
818 | //
|
---|
819 | Int_t MPointingDevCalc::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
|
---|
820 | {
|
---|
821 | Bool_t rc = kFALSE;
|
---|
822 | if (IsEnvDefined(env, prefix, "NumMinStars", print))
|
---|
823 | {
|
---|
824 | SetNumMinStars(GetEnvValue(env, prefix, "NumMinStars", (Int_t)fNumMinStars));
|
---|
825 | rc = kTRUE;
|
---|
826 | }
|
---|
827 | if (IsEnvDefined(env, prefix, "NsbLevel", print))
|
---|
828 | {
|
---|
829 | SetNsbLevel(GetEnvValue(env, prefix, "NsbLevel", fNsbLevel));
|
---|
830 | rc = kTRUE;
|
---|
831 | }
|
---|
832 | if (IsEnvDefined(env, prefix, "NsbMin", print))
|
---|
833 | {
|
---|
834 | SetNsbMin(GetEnvValue(env, prefix, "NsbMin", fNsbMin));
|
---|
835 | rc = kTRUE;
|
---|
836 | }
|
---|
837 | if (IsEnvDefined(env, prefix, "NsbMax", print))
|
---|
838 | {
|
---|
839 | SetNsbMax(GetEnvValue(env, prefix, "NsbMax", fNsbMax));
|
---|
840 | rc = kTRUE;
|
---|
841 | }
|
---|
842 | if (IsEnvDefined(env, prefix, "MaxAbsDev", print))
|
---|
843 | {
|
---|
844 | SetMaxAbsDev(GetEnvValue(env, prefix, "MaxAbsDev", fMaxAbsDev));
|
---|
845 | rc = kTRUE;
|
---|
846 | }
|
---|
847 | if (IsEnvDefined(env, prefix, "Dx", print))
|
---|
848 | {
|
---|
849 | fDx = GetEnvValue(env, prefix, "Dx", fDx);
|
---|
850 | rc = kTRUE;
|
---|
851 | }
|
---|
852 | if (IsEnvDefined(env, prefix, "Dy", print))
|
---|
853 | {
|
---|
854 | fDy = GetEnvValue(env, prefix, "Dy", fDy);
|
---|
855 | rc = kTRUE;
|
---|
856 | }
|
---|
857 | if (IsEnvDefined(env, prefix, "MaxAge", print))
|
---|
858 | {
|
---|
859 | fMaxAge = GetEnvValue(env, prefix, "MaxAge", fMaxAge);
|
---|
860 | rc = kTRUE;
|
---|
861 | }
|
---|
862 | if (IsEnvDefined(env, prefix, "FilePrefix", print))
|
---|
863 | {
|
---|
864 | fFilePrefix = GetEnvValue(env, prefix, "FilePrefix", fFilePrefix);
|
---|
865 | rc = kTRUE;
|
---|
866 | }
|
---|
867 | if (IsEnvDefined(env, prefix, "PointingModels", print))
|
---|
868 | {
|
---|
869 | SetPointingModels(GetEnvValue(env, prefix, "PointingModels", GetPointingModels()));
|
---|
870 | rc = kTRUE;
|
---|
871 | }
|
---|
872 |
|
---|
873 | return rc;
|
---|
874 | }
|
---|