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:
|
---|
119 | //
|
---|
120 | // 18.03.2006: The camera holding has been repaired and the camera got
|
---|
121 | // shifted a little bit.
|
---|
122 | //
|
---|
123 | // 16.04.2006: Around this date a roque lamp focussing hass been done
|
---|
124 | //
|
---|
125 | // 25.04.2006: A missalignment intrduced with the roque adjust has been
|
---|
126 | // corrected
|
---|
127 | //
|
---|
128 | // The starguider pointing model for the time before 18.3.2006 and after
|
---|
129 | // April 2006 (in fact there are no TPoints until 07/2006 to check it)
|
---|
130 | // and for the period 07/2006 to (at least) 06/2007 are very similar.
|
---|
131 | //
|
---|
132 | // The pointing model for the time between 18.3.2006 and 04/2006 is
|
---|
133 | // clearly different, mainly giving different Azimuth values between
|
---|
134 | // Zenith and roughly ~25deg, and a slight offset on both axes.
|
---|
135 | //
|
---|
136 | // 10.5.2006: pos1 -= pos0 commented (What was the mentioned fix?)
|
---|
137 | // 29.6.2006: repaired
|
---|
138 | //
|
---|
139 | // 23.3.2006: new starguider algorithm
|
---|
140 | //
|
---|
141 | // 17.3.2005: Fixed units of "nompos" in MDriveCom
|
---|
142 | //
|
---|
143 | //
|
---|
144 | // New pointing models have been installed (if the pointing model
|
---|
145 | // is different, than the previous one it might mean, that also
|
---|
146 | // the starguider calibration is different.) The date only means
|
---|
147 | // day-time (noon) at which the model has been changed.
|
---|
148 | //
|
---|
149 | // 29. Apr. 2004 ~25800
|
---|
150 | // 5. Aug. 2004 ~32000
|
---|
151 | // 19. Aug. 2004 ~33530
|
---|
152 | // 7. Jun. 2005 ~57650
|
---|
153 | // 8. Jun. 2005
|
---|
154 | // 9. Jun. 2005 ~57860
|
---|
155 | // 12. Sep. 2005 ~68338
|
---|
156 | // 24. Nov. 2005 ~75562
|
---|
157 | // 17. Oct. 2006 ~103130
|
---|
158 | // 17. Jun. 2007 ~248193
|
---|
159 | // 18. Oct. 2007 291104 // Correction for the offsets introduced by AMC
|
---|
160 | // 14. Jan. 2008 328198 // Complete new pointing model
|
---|
161 | // 11. Jun. 2008 1000534 (ca. 23:00) // Before new TPoints
|
---|
162 | // 19. Jun. 2008 (ca. 15:00) // From new TPoints
|
---|
163 | // 7. Mar. 2009 (ca. 14:00) // From new TPoints (0808-0902)
|
---|
164 | //
|
---|
165 | // From 2.2.2006 beginnig of the night (21:05h, run >=81855) to 24.2.2006
|
---|
166 | // beginning of the the night (20:34h, run<83722) the LEDs did not work.
|
---|
167 | // In the time after this incident the shift crew often forgot to switch on
|
---|
168 | // the LEDs at the beginning of the night!
|
---|
169 | //
|
---|
170 | // [2006-03-07 00:10:58] In the daytime, we raised the position of the
|
---|
171 | // 9 o'clock LED by one screw hole to make it visible when the TPoint
|
---|
172 | // Lid is closed. (< run 84980)
|
---|
173 | //
|
---|
174 | // Mirror refocussing around 23.Apr.2006, from the Runbook.
|
---|
175 | //
|
---|
176 | // 25./26.4.2006: Run 89180
|
---|
177 | //
|
---|
178 | // (Note: The year here is not a typo!)
|
---|
179 | // [2007-04-25 23:50:39]
|
---|
180 | // Markus is performing AMC focussing.
|
---|
181 | //
|
---|
182 | // Mirror refocussing around 4.Aug.2007, from the Runbook:
|
---|
183 | //
|
---|
184 | // [2007-08-04 04:46:47]
|
---|
185 | // We finished with the focussing with Polaris. The images need to be
|
---|
186 | // analysed and new LUTs generated.
|
---|
187 | //
|
---|
188 | // [2007-08-04 23:47:30]
|
---|
189 | // Actually we see that the mispointing is always large; probably since
|
---|
190 | // the LUT tables have not yet been adjusted to the new focussing.
|
---|
191 | //
|
---|
192 | // [2007-08-03 23:07:58]
|
---|
193 | // Data taking stopped. Mirror focussing.
|
---|
194 | //
|
---|
195 | // [2007-08-05 00:09:16]
|
---|
196 | // We take some pictures on stars nearby Cyg X3 with the sbig camera;
|
---|
197 | // actually the spot doesn't look very nice... The pictures have been
|
---|
198 | // saved with name Deneb- and Sadr- Polaris seems a bit better. Should we
|
---|
199 | // have new LUT tables after the focussing?
|
---|
200 | //
|
---|
201 | // [2007-08-10 20:18:48]
|
---|
202 | // Tonight we take first images of Polaris with a new LUT file generated
|
---|
203 | // based on the recent focussing. The image will be analysed tomorrow and
|
---|
204 | // than new LUTs will be generated. For tonight the focussing is still not
|
---|
205 | // changed.
|
---|
206 | //
|
---|
207 | // [2007-08-14 20:57:59]
|
---|
208 | // The weather is fine. There is a group of hobby astronomers at the
|
---|
209 | // helicopter parking. Before data taking, we tried to check the new LUTs.
|
---|
210 | // However, because of technical problems with the new LUTs we had to
|
---|
211 | // postpone the measurements. We lost some 10 min of data taking because
|
---|
212 | // of this.
|
---|
213 | //
|
---|
214 | // [2007-08-14 22:29:37] Run 267253
|
---|
215 | // Before continuing the observation we perform a focussing test with the
|
---|
216 | // new LUTs from recent Polaris focussing. Note: Data on Her X-1 was taken
|
---|
217 | // with old focussing. We performed PSF measurements on Kornephorus (Zd
|
---|
218 | // 31.77, Az 264,67) first with old LUTs and then with new LUTs. We took
|
---|
219 | // T-points with both focussing. The first T-Point corresponds to the
|
---|
220 | // previous focussing and the second with the improved. Please check both
|
---|
221 | // T-points for eventual misspointing. We found big improvement of the PSF
|
---|
222 | // with the new LUTs and will therefore continue from now on with the new
|
---|
223 | // focussing. Having a first look at the SBIG pictures we see a slightly
|
---|
224 | // misspointing of ~0.1 deg with the new LUTs.
|
---|
225 | //
|
---|
226 | // [2007-08-14 22:46:10]
|
---|
227 | // Comparing the trigger rate with yesterday night we do not see an
|
---|
228 | // improvement with the new focussing.
|
---|
229 | //
|
---|
230 | // [2007-10-27 email]
|
---|
231 | // [...]
|
---|
232 | // When removing the cover of the motor all 4 nuts which fix the gear of
|
---|
233 | // the motor to the structure fell out. The motor was completely loose.
|
---|
234 | // [...]
|
---|
235 | // We checked also the second azimuth motor and .. the same situation.
|
---|
236 | // [...]
|
---|
237 | // Peter Sawallisch opened and fixed all 3 motors of MAGIC-I. We checked
|
---|
238 | // all critical screws that came to our minds and fixed them as much as
|
---|
239 | // possible with loctite and counter nuts. No damage whatsoever has been
|
---|
240 | // found.
|
---|
241 | //
|
---|
242 | //
|
---|
243 | // ToDo:
|
---|
244 | // -----
|
---|
245 | //
|
---|
246 | // * Is 0/0 the best assumption if the starguider partly fails?
|
---|
247 | //
|
---|
248 | //
|
---|
249 | // Input Container:
|
---|
250 | // MRawRunHeader
|
---|
251 | // MReportStarguider
|
---|
252 | //
|
---|
253 | // Output Container:
|
---|
254 | // MPointingDev
|
---|
255 | //
|
---|
256 | /////////////////////////////////////////////////////////////////////////////
|
---|
257 | #include "MPointingDevCalc.h"
|
---|
258 |
|
---|
259 | #include <stdlib.h> // atoi (Ubuntu 8.10)
|
---|
260 |
|
---|
261 | #include "MLog.h"
|
---|
262 | #include "MLogManip.h"
|
---|
263 |
|
---|
264 | #include "MParList.h"
|
---|
265 |
|
---|
266 | #include "MAstro.h"
|
---|
267 | #include "MPointing.h"
|
---|
268 | #include "MPointingDev.h"
|
---|
269 | #include "MRawRunHeader.h"
|
---|
270 | #include "MReportStarguider.h"
|
---|
271 |
|
---|
272 | ClassImp(MPointingDevCalc);
|
---|
273 |
|
---|
274 | using namespace std;
|
---|
275 |
|
---|
276 | const TString MPointingDevCalc::fgFilePrefix="resources/starguider";
|
---|
277 |
|
---|
278 | // --------------------------------------------------------------------------
|
---|
279 | //
|
---|
280 | // Destructor. Call Clear() and delete fPointingModels if any.
|
---|
281 | //
|
---|
282 | MPointingDevCalc::~MPointingDevCalc()
|
---|
283 | {
|
---|
284 | Clear();
|
---|
285 | }
|
---|
286 |
|
---|
287 | // --------------------------------------------------------------------------
|
---|
288 | //
|
---|
289 | // Delete fPointing and set pointing to NULL
|
---|
290 | //
|
---|
291 | void MPointingDevCalc::Clear(Option_t *o)
|
---|
292 | {
|
---|
293 | if (fPointing)
|
---|
294 | delete fPointing;
|
---|
295 |
|
---|
296 | fPointing = NULL;
|
---|
297 | }
|
---|
298 |
|
---|
299 | // --------------------------------------------------------------------------
|
---|
300 | //
|
---|
301 | // Sort the entries in fPoinitngModels
|
---|
302 | //
|
---|
303 | void MPointingDevCalc::SortPointingModels()
|
---|
304 | {
|
---|
305 | const int n = fPointingModels.GetSize();
|
---|
306 |
|
---|
307 | TArrayI idx(n);
|
---|
308 |
|
---|
309 | TMath::Sort(n, fPointingModels.GetArray(), idx.GetArray(), kFALSE);
|
---|
310 |
|
---|
311 | const TArrayI arr(fPointingModels);
|
---|
312 |
|
---|
313 | for (int i=0; i<n; i++)
|
---|
314 | fPointingModels[i] = arr[idx[i]];
|
---|
315 | }
|
---|
316 |
|
---|
317 | // --------------------------------------------------------------------------
|
---|
318 | //
|
---|
319 | // Set new pointing models
|
---|
320 | //
|
---|
321 | void MPointingDevCalc::SetPointingModels(const TString &models)
|
---|
322 | {
|
---|
323 | fPointingModels.Set(0);
|
---|
324 |
|
---|
325 | if (models.IsNull())
|
---|
326 | return;
|
---|
327 |
|
---|
328 | TObjArray *arr = models.Tokenize(" ");
|
---|
329 |
|
---|
330 | const int n = arr->GetEntries();
|
---|
331 | fPointingModels.Set(n);
|
---|
332 |
|
---|
333 | for (int i=0; i<n; i++)
|
---|
334 | fPointingModels[i] = atoi((*arr)[i]->GetName());
|
---|
335 |
|
---|
336 | delete arr;
|
---|
337 |
|
---|
338 | SortPointingModels();
|
---|
339 | }
|
---|
340 |
|
---|
341 | // --------------------------------------------------------------------------
|
---|
342 | //
|
---|
343 | // Return a string with the pointing models, seperated by a space.
|
---|
344 | //
|
---|
345 | TString MPointingDevCalc::GetPointingModels() const
|
---|
346 | {
|
---|
347 | TString rc;
|
---|
348 | for (int i=0; i<fPointingModels.GetSize(); i++)
|
---|
349 | rc += Form ("%d ", fPointingModels[i]);
|
---|
350 |
|
---|
351 | return rc;
|
---|
352 | }
|
---|
353 |
|
---|
354 | // --------------------------------------------------------------------------
|
---|
355 | //
|
---|
356 | // Add a number to the pointing models
|
---|
357 | //
|
---|
358 | void MPointingDevCalc::AddPointingModel(UInt_t runnum)
|
---|
359 | {
|
---|
360 | const int n = fPointingModels.GetSize();
|
---|
361 | for (int i=0; i<n; i++)
|
---|
362 | if ((UInt_t)fPointingModels[i]==runnum)
|
---|
363 | {
|
---|
364 | *fLog << warn << "WARNING - Pointing model " << runnum << " already in list... ignored." << endl;
|
---|
365 | return;
|
---|
366 | }
|
---|
367 |
|
---|
368 | fPointingModels.Set(n+1);
|
---|
369 | fPointingModels[n] = runnum;
|
---|
370 |
|
---|
371 | SortPointingModels();
|
---|
372 | }
|
---|
373 |
|
---|
374 | // --------------------------------------------------------------------------
|
---|
375 | //
|
---|
376 | // Find the highest number in the array which is lower or equal num.
|
---|
377 | //
|
---|
378 | UInt_t MPointingDevCalc::FindPointingModel(UInt_t num)
|
---|
379 | {
|
---|
380 | const int n = fPointingModels.GetSize();
|
---|
381 | if (n==0)
|
---|
382 | return (UInt_t)-1;
|
---|
383 |
|
---|
384 | // Loop over all pointing models
|
---|
385 | for (int i=0; i<n; i++)
|
---|
386 | {
|
---|
387 | // The number stored in the array did not yet overtake the runnumber
|
---|
388 | if ((UInt_t)fPointingModels[i]<num)
|
---|
389 | continue;
|
---|
390 |
|
---|
391 | // The first pointing model is later than this run: use a default
|
---|
392 | if (i==0)
|
---|
393 | return (UInt_t)-1;
|
---|
394 |
|
---|
395 | // The last entry in the array is the right one.
|
---|
396 | return fPointingModels[i-1];
|
---|
397 | }
|
---|
398 |
|
---|
399 | // Runnumber is after last entry of pointing models. Use the last one.
|
---|
400 | return fPointingModels[n-1];
|
---|
401 | }
|
---|
402 |
|
---|
403 | // --------------------------------------------------------------------------
|
---|
404 | //
|
---|
405 | // Clear the pointing model. If run-number >= 87751 read the new
|
---|
406 | // pointing model with fFilePrefix
|
---|
407 | //
|
---|
408 | Bool_t MPointingDevCalc::ReadPointingModel(const MRawRunHeader &run)
|
---|
409 | {
|
---|
410 | const UInt_t num = FindPointingModel(run.GetRunNumber());
|
---|
411 |
|
---|
412 | // No poinitng models are defined. Use simple dx/dy-calibration
|
---|
413 | if (num==(UInt_t)-1)
|
---|
414 | {
|
---|
415 | Clear();
|
---|
416 | return kTRUE;
|
---|
417 | }
|
---|
418 |
|
---|
419 | // compile the name for the starguider files
|
---|
420 | // The file with the number 00000000 is the default file
|
---|
421 | TString fname = Form("%s%08d.txt", fFilePrefix.Data(), num);
|
---|
422 |
|
---|
423 | if (!fPointing)
|
---|
424 | fPointing = new MPointing;
|
---|
425 |
|
---|
426 | if (fname==fPointing->GetName())
|
---|
427 | {
|
---|
428 | *fLog << inf << fname << " already loaded." << endl;
|
---|
429 | return kTRUE;
|
---|
430 | }
|
---|
431 |
|
---|
432 | return fPointing->Load(fname);
|
---|
433 | }
|
---|
434 |
|
---|
435 | // --------------------------------------------------------------------------
|
---|
436 | //
|
---|
437 | // Check the file/run type from the run-header and if it is a data file
|
---|
438 | // load starguider calibration.
|
---|
439 | //
|
---|
440 | Bool_t MPointingDevCalc::ReInit(MParList *plist)
|
---|
441 | {
|
---|
442 | MRawRunHeader *run = (MRawRunHeader*)plist->FindObject("MRawRunHeader");
|
---|
443 | if (!run)
|
---|
444 | {
|
---|
445 | *fLog << err << "MRawRunHeader not found... aborting." << endl;
|
---|
446 | return kFALSE;
|
---|
447 | }
|
---|
448 |
|
---|
449 | fNsbSum = 0;
|
---|
450 | fNsbSq = 0;
|
---|
451 | fNsbCount = 0;
|
---|
452 |
|
---|
453 | fRunType = run->GetRunType();
|
---|
454 |
|
---|
455 | switch (fRunType)
|
---|
456 | {
|
---|
457 | case MRawRunHeader::kRTData:
|
---|
458 | if (!fReport)
|
---|
459 | *fLog << warn << "MReportStarguider not found... skipped." << endl;
|
---|
460 | return ReadPointingModel(*run);
|
---|
461 |
|
---|
462 | case MRawRunHeader::kRTMonteCarlo:
|
---|
463 | return kTRUE;
|
---|
464 |
|
---|
465 | case MRawRunHeader::kRTPedestal:
|
---|
466 | *fLog << err << "Cannot work in a pedestal Run!... aborting." << endl;
|
---|
467 | return kFALSE;
|
---|
468 |
|
---|
469 | case MRawRunHeader::kRTCalibration:
|
---|
470 | *fLog << err << "Cannot work in a calibration Run!... aborting." << endl;
|
---|
471 | return kFALSE;
|
---|
472 |
|
---|
473 | default:
|
---|
474 | *fLog << err << "Run Type " << fRunType << " unknown!... aborting." << endl;
|
---|
475 | return kFALSE;
|
---|
476 | }
|
---|
477 |
|
---|
478 | return kTRUE;
|
---|
479 | }
|
---|
480 |
|
---|
481 | // --------------------------------------------------------------------------
|
---|
482 | //
|
---|
483 | // Search for 'MPointingPos'. Create if not found.
|
---|
484 | //
|
---|
485 | Int_t MPointingDevCalc::PreProcess(MParList *plist)
|
---|
486 | {
|
---|
487 | fDeviation = (MPointingDev*)plist->FindCreateObj("MPointingDev");
|
---|
488 | fReport = (MReportStarguider*)plist->FindObject("MReportStarguider");
|
---|
489 |
|
---|
490 | // We use kRTNone here as a placeholder for data runs.
|
---|
491 | fRunType = MRawRunHeader::kRTNone;
|
---|
492 | fLastMjd = -1;
|
---|
493 |
|
---|
494 | fSkip.Reset();
|
---|
495 |
|
---|
496 | // In cases in which ReInit isn't called right in time (e.g. if
|
---|
497 | // the first starguider event comes before the first data event)
|
---|
498 | fNsbSum = 0;
|
---|
499 | fNsbSq = 0;
|
---|
500 | fNsbCount = 0;
|
---|
501 |
|
---|
502 | return fDeviation ? kTRUE : kFALSE;
|
---|
503 | }
|
---|
504 |
|
---|
505 | // --------------------------------------------------------------------------
|
---|
506 | //
|
---|
507 | // Increase fSkip[i] by one. If the data in fDeviation is outdated (older
|
---|
508 | // than fMaxAge) and the current report should be skipped reset DevZdAz and
|
---|
509 | // DevXY and fSkip[6] is increased by one.
|
---|
510 | //
|
---|
511 | void MPointingDevCalc::Skip(Int_t i)
|
---|
512 | {
|
---|
513 | fSkip[i]++;
|
---|
514 |
|
---|
515 | const Double_t diff = (fReport->GetMjd()-fLastMjd)*1440; // [min] 1440=24*60
|
---|
516 | if (diff<fMaxAge && fLastMjd>0)
|
---|
517 | return;
|
---|
518 |
|
---|
519 | fDeviation->SetDevZdAz(0, 0);
|
---|
520 | fDeviation->SetDevXY(0, 0);
|
---|
521 | fSkip[6]++;
|
---|
522 | }
|
---|
523 |
|
---|
524 | // --------------------------------------------------------------------------
|
---|
525 | //
|
---|
526 | // Do a full starguider calibration using a pointing model for the starguider.
|
---|
527 | //
|
---|
528 | void MPointingDevCalc::DoCalibration(Double_t devzd, Double_t devaz) const
|
---|
529 | {
|
---|
530 | if (!fPointing)
|
---|
531 | {
|
---|
532 | // Do a simple starguider calibration using a simple offset in x and y
|
---|
533 | fDeviation->SetDevZdAz(devzd, devaz);
|
---|
534 |
|
---|
535 | // Linear starguider calibration taken from April/May data
|
---|
536 | // For calibration add MDriveReport::GetErrorZd/Az !
|
---|
537 | fDeviation->SetDevXY(fDx, fDy); // 1arcmin ~ 5mm
|
---|
538 |
|
---|
539 | return;
|
---|
540 | }
|
---|
541 |
|
---|
542 | // Get the nominal position the star is at the sky
|
---|
543 | // Unit: deg
|
---|
544 | ZdAz nom(fReport->GetNominalZd(), fReport->GetNominalAz());
|
---|
545 | nom *= TMath::DegToRad();
|
---|
546 |
|
---|
547 | // Get the mispointing measured by the telescope. It is
|
---|
548 | // calculate as follows:
|
---|
549 | //
|
---|
550 | // The mispointing measured by the starguider:
|
---|
551 | // ZdAz mis(devzd, devaz);
|
---|
552 | // mis *= TMath::DegToRad();
|
---|
553 |
|
---|
554 | // The pointing model is the conversion from the real pointing
|
---|
555 | // position of the telescope into the pointing position measured
|
---|
556 | // by the starguider.
|
---|
557 | //
|
---|
558 | // To keep as close to the fitted model we use the forward correction.
|
---|
559 |
|
---|
560 | // Position at which the starguider camera is pointing in real:
|
---|
561 | // pointing position = nominal position - dev
|
---|
562 | //
|
---|
563 | // The position measured as the starguider's pointing position
|
---|
564 | ZdAz pos(nom); // cpos = sao - dev
|
---|
565 | pos -= ZdAz(devzd, devaz)*TMath::DegToRad();
|
---|
566 |
|
---|
567 | // Now we convert the starguider's pointing position into the
|
---|
568 | // telescope pointing position (the pointing model converts
|
---|
569 | // the telescope pointing position into the starguider pointing
|
---|
570 | // position)
|
---|
571 | ZdAz point = fPointing->CorrectBack(pos); //FWD!!!
|
---|
572 |
|
---|
573 | // MSrcPosCalc uses the following condition to calculate the
|
---|
574 | // source position in the camera:
|
---|
575 | // real pointing pos = nominal pointing pos - dev
|
---|
576 | // --> dev = nominal - real
|
---|
577 | // Therefor we calculate dev as follows:
|
---|
578 | ZdAz dev(nom);
|
---|
579 | dev -= point;
|
---|
580 | dev *= TMath::RadToDeg();
|
---|
581 |
|
---|
582 | /*
|
---|
583 | // We chose the other way. It is less accurate because is is the
|
---|
584 | // other was than the poinitng model was fittet, but it is more
|
---|
585 | // accurate because the nominal (i.e. real) pointing position
|
---|
586 | // is less accurately known than the position returned by the
|
---|
587 | // starguider.
|
---|
588 | //
|
---|
589 | // Calculate the deviation which would be measured by the starguider
|
---|
590 | // if applied to a perfectly pointing telescope.
|
---|
591 | ZdAz dev = fPointing->Correct(nom);
|
---|
592 | dev -= nom;
|
---|
593 |
|
---|
594 | // Now add these offsets and the starguider measured offsets to
|
---|
595 | // the real pointing deviation of the telescope (note, that
|
---|
596 | // signs here are just conventions)
|
---|
597 | dev += ZdAz(devzd, devaz)*TMath::DegToRad(); // --> nom-mis
|
---|
598 | dev *= TMath::RadToDeg();
|
---|
599 | */
|
---|
600 |
|
---|
601 | // Check if the starguider pointing model requests overwriting
|
---|
602 | // of the values with constants (e.g. 0)
|
---|
603 | devaz = fPointing->IsPxValid() ? fPointing->GetPx() : dev.Az();
|
---|
604 | devzd = fPointing->IsPyValid() ? fPointing->GetPy() : dev.Zd();
|
---|
605 |
|
---|
606 | fDeviation->SetDevZdAz(devzd, devaz);
|
---|
607 | fDeviation->SetDevXY(fPointing->GetDxy());
|
---|
608 | }
|
---|
609 |
|
---|
610 | Int_t MPointingDevCalc::ProcessStarguiderReport()
|
---|
611 | {
|
---|
612 | Double_t devzd = fReport->GetDevZd(); // [arcmin]
|
---|
613 | Double_t devaz = fReport->GetDevAz(); // [arcmin]
|
---|
614 | if (devzd==0 && devaz==0)
|
---|
615 | {
|
---|
616 | Skip(1);
|
---|
617 | return kTRUE;
|
---|
618 | }
|
---|
619 |
|
---|
620 | if (!fReport->IsMonitoring())
|
---|
621 | {
|
---|
622 | Skip(2);
|
---|
623 | return kTRUE;
|
---|
624 | }
|
---|
625 |
|
---|
626 | devzd /= 60; // Convert from arcmin to deg
|
---|
627 | devaz /= 60; // Convert from arcmin to deg
|
---|
628 |
|
---|
629 | const Double_t nsb = fReport->GetSkyBrightness();
|
---|
630 | if (nsb>0)
|
---|
631 | {
|
---|
632 | if (nsb>fNsbMin && nsb<fNsbMax)
|
---|
633 | {
|
---|
634 | fNsbSum += nsb;
|
---|
635 | fNsbSq += nsb*nsb;
|
---|
636 | fNsbCount++;
|
---|
637 | }
|
---|
638 |
|
---|
639 | if (fNsbCount>0)
|
---|
640 | {
|
---|
641 | const Double_t sum = fNsbSum/fNsbCount;
|
---|
642 | const Double_t sq = fNsbSq /fNsbCount;
|
---|
643 |
|
---|
644 | const Double_t rms = fNsbLevel*TMath::Sqrt(sq - sum*sum);
|
---|
645 |
|
---|
646 | if (nsb<sum-rms || nsb>sum+rms)
|
---|
647 | {
|
---|
648 | Skip(3);
|
---|
649 | return kTRUE;
|
---|
650 | }
|
---|
651 | }
|
---|
652 |
|
---|
653 | if (fReport->GetNumIdentifiedStars()<fNumMinStars)
|
---|
654 | {
|
---|
655 | Skip(4);
|
---|
656 | return kTRUE;
|
---|
657 | }
|
---|
658 | }
|
---|
659 |
|
---|
660 | // >= 87751 (31.3.06 12:00)
|
---|
661 |
|
---|
662 | // Calculate absolute deviation
|
---|
663 | const Double_t dev = MAstro::GetDevAbs(fReport->GetNominalZd(), devzd, devaz);
|
---|
664 |
|
---|
665 | // Sanity check... larger deviation are strange and ignored
|
---|
666 | if (dev*60>fMaxAbsDev)
|
---|
667 | {
|
---|
668 | Skip(5);
|
---|
669 | return kTRUE;
|
---|
670 | }
|
---|
671 |
|
---|
672 | DoCalibration(devzd, devaz);
|
---|
673 |
|
---|
674 | fSkip[0]++;
|
---|
675 | fLastMjd = fReport->GetMjd();
|
---|
676 |
|
---|
677 | return kTRUE;
|
---|
678 | }
|
---|
679 |
|
---|
680 | // --------------------------------------------------------------------------
|
---|
681 | //
|
---|
682 | // See class description.
|
---|
683 | //
|
---|
684 | Int_t MPointingDevCalc::Process()
|
---|
685 | {
|
---|
686 | switch (fRunType)
|
---|
687 | {
|
---|
688 | case MRawRunHeader::kRTNone:
|
---|
689 | case MRawRunHeader::kRTData:
|
---|
690 | return fReport ? ProcessStarguiderReport() : kTRUE;
|
---|
691 |
|
---|
692 | case MRawRunHeader::kRTMonteCarlo:
|
---|
693 | fSkip[0]++;
|
---|
694 | fDeviation->SetDevZdAz(0, 0);
|
---|
695 | fDeviation->SetDevXY(0, 0);
|
---|
696 | return kTRUE;
|
---|
697 | }
|
---|
698 | return kTRUE;
|
---|
699 | }
|
---|
700 |
|
---|
701 | // --------------------------------------------------------------------------
|
---|
702 | //
|
---|
703 | // Print execution statistics
|
---|
704 | //
|
---|
705 | Int_t MPointingDevCalc::PostProcess()
|
---|
706 | {
|
---|
707 | if (GetNumExecutions()==0)
|
---|
708 | return kTRUE;
|
---|
709 |
|
---|
710 | *fLog << inf << endl;
|
---|
711 | *fLog << GetDescriptor() << " execution statistics:" << endl;
|
---|
712 | PrintSkipped(fSkip[1], "Starguider deviation not set, is exactly 0/0");
|
---|
713 | PrintSkipped(fSkip[2], "Starguider was not monitoring (eg. LEDs off)");
|
---|
714 | PrintSkipped(fSkip[3], Form("NSB out of %.1f sigma range", fNsbLevel));
|
---|
715 | PrintSkipped(fSkip[4], Form("Number of identified stars < %d", fNumMinStars));
|
---|
716 | PrintSkipped(fSkip[5], Form("Absolute deviation > %.1farcmin", fMaxAbsDev));
|
---|
717 | PrintSkipped(fSkip[6], Form("Events set to 0 because older than %.1fmin", fMaxAge));
|
---|
718 | *fLog << " " << (int)fSkip[0] << " (" << Form("%5.1f", 100.*fSkip[0]/GetNumExecutions()) << "%) Evts survived calculation!" << endl;
|
---|
719 | *fLog << endl;
|
---|
720 |
|
---|
721 | return kTRUE;
|
---|
722 | }
|
---|
723 |
|
---|
724 | // --------------------------------------------------------------------------
|
---|
725 | //
|
---|
726 | // MPointingDevCalc.NumMinStars: 8
|
---|
727 | // MPointingDevCalc.NsbLevel: 3.0
|
---|
728 | // MPointingDevCalc.NsbMin: 30
|
---|
729 | // MPointingDevCalc.NsbMax: 60
|
---|
730 | // MPointingDevCalc.MaxAbsDev: 15
|
---|
731 | // MPointingDevCalc.MaxAge: 1.0
|
---|
732 | // MPointingDevCalc.Dx: -0.001
|
---|
733 | // MPointingDevCalc.Dy: -0.004
|
---|
734 | //
|
---|
735 | // For a detailed description see the class reference.
|
---|
736 | //
|
---|
737 | Int_t MPointingDevCalc::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
|
---|
738 | {
|
---|
739 | Bool_t rc = kFALSE;
|
---|
740 | if (IsEnvDefined(env, prefix, "NumMinStars", print))
|
---|
741 | {
|
---|
742 | SetNumMinStars(GetEnvValue(env, prefix, "NumMinStars", (Int_t)fNumMinStars));
|
---|
743 | rc = kTRUE;
|
---|
744 | }
|
---|
745 | if (IsEnvDefined(env, prefix, "NsbLevel", print))
|
---|
746 | {
|
---|
747 | SetNsbLevel(GetEnvValue(env, prefix, "NsbLevel", fNsbLevel));
|
---|
748 | rc = kTRUE;
|
---|
749 | }
|
---|
750 | if (IsEnvDefined(env, prefix, "NsbMin", print))
|
---|
751 | {
|
---|
752 | SetNsbMin(GetEnvValue(env, prefix, "NsbMin", fNsbMin));
|
---|
753 | rc = kTRUE;
|
---|
754 | }
|
---|
755 | if (IsEnvDefined(env, prefix, "NsbMax", print))
|
---|
756 | {
|
---|
757 | SetNsbMax(GetEnvValue(env, prefix, "NsbMax", fNsbMax));
|
---|
758 | rc = kTRUE;
|
---|
759 | }
|
---|
760 | if (IsEnvDefined(env, prefix, "MaxAbsDev", print))
|
---|
761 | {
|
---|
762 | SetMaxAbsDev(GetEnvValue(env, prefix, "MaxAbsDev", fMaxAbsDev));
|
---|
763 | rc = kTRUE;
|
---|
764 | }
|
---|
765 | if (IsEnvDefined(env, prefix, "Dx", print))
|
---|
766 | {
|
---|
767 | fDx = GetEnvValue(env, prefix, "Dx", fDx);
|
---|
768 | rc = kTRUE;
|
---|
769 | }
|
---|
770 | if (IsEnvDefined(env, prefix, "Dy", print))
|
---|
771 | {
|
---|
772 | fDy = GetEnvValue(env, prefix, "Dy", fDy);
|
---|
773 | rc = kTRUE;
|
---|
774 | }
|
---|
775 | if (IsEnvDefined(env, prefix, "MaxAge", print))
|
---|
776 | {
|
---|
777 | fMaxAge = GetEnvValue(env, prefix, "MaxAge", fMaxAge);
|
---|
778 | rc = kTRUE;
|
---|
779 | }
|
---|
780 | if (IsEnvDefined(env, prefix, "FilePrefix", print))
|
---|
781 | {
|
---|
782 | fFilePrefix = GetEnvValue(env, prefix, "FilePrefix", fFilePrefix);
|
---|
783 | rc = kTRUE;
|
---|
784 | }
|
---|
785 | if (IsEnvDefined(env, prefix, "PointingModels", print))
|
---|
786 | {
|
---|
787 | SetPointingModels(GetEnvValue(env, prefix, "PointingModels", GetPointingModels()));
|
---|
788 | rc = kTRUE;
|
---|
789 | }
|
---|
790 |
|
---|
791 | return rc;
|
---|
792 | }
|
---|