source: trunk/MagicSoft/Cosy/tpoint/TPointStar.cc@ 9029

Last change on this file since 9029 was 8823, checked in by tbretz, 17 years ago
*** empty log message ***
File size: 4.6 KB
Line 
1#include "TPointStar.h"
2
3#include <TMath.h>
4
5#include "MLog.h"
6#include "MLogManip.h"
7
8#include "MPointing.h"
9
10ClassImp(TPointStar);
11
12using namespace std;
13
14TPointStar::TPointStar(Double_t sel, Double_t saz, Double_t rel, Double_t raz) :
15 fStarAz(saz*TMath::DegToRad()),
16 fStarEl(sel*TMath::DegToRad()),
17 fRawAz(raz*TMath::DegToRad()),
18 fRawEl(rel*TMath::DegToRad()), fMag(-25)
19{
20}
21
22Double_t TPointStar::GetDEl() const { return (fRawEl-fStarEl)*TMath::RadToDeg(); }
23Double_t TPointStar::GetDZd() const { return -GetDEl(); }
24Double_t TPointStar::GetDAz() const { return (fRawAz-fStarAz)*TMath::RadToDeg(); }
25Double_t TPointStar::GetStarEl() const { return fStarEl*TMath::RadToDeg(); }
26Double_t TPointStar::GetStarZd() const { return 90.-fStarEl*TMath::RadToDeg(); }
27Double_t TPointStar::GetStarAz() const { return fStarAz*TMath::RadToDeg(); }
28Double_t TPointStar::GetRawEl() const { return fRawEl*TMath::RadToDeg(); }
29Double_t TPointStar::GetRawAz() const { return fRawAz*TMath::RadToDeg(); }
30Double_t TPointStar::GetRawZd() const { return 90.-fRawEl*TMath::RadToDeg(); }
31
32ZdAz TPointStar::GetStarZdAz() const { return ZdAz(TMath::Pi()/2-fStarEl, fStarAz); }
33AltAz TPointStar::GetStarAltAz() const { return AltAz(fStarEl, fStarAz); }
34
35ZdAz TPointStar::GetRawZdAz() const { return ZdAz(TMath::Pi()/2-fRawEl, fRawAz); }
36AltAz TPointStar::GetRawAltAz() const { return AltAz(fRawEl, fRawAz); }
37
38void TPointStar::AdjustEl(Double_t del) { fStarEl += del*TMath::DegToRad(); }
39void TPointStar::AdjustAz(Double_t daz) { fStarAz += daz*TMath::DegToRad(); }
40
41void TPointStar::Adjust(const MPointing &bend)
42{
43 AltAz p = bend(GetStarAltAz());
44 fStarEl = p.Alt();
45 fStarAz = p.Az();
46}
47
48void TPointStar::AdjustBack(const MPointing &bend)
49{
50 AltAz p = bend.CorrectBack(GetRawAltAz());
51 fRawEl = p.Alt();
52 fRawAz = p.Az();
53}
54
55Double_t TPointStar::GetResidual(Double_t *err) const
56 {
57 /*
58 TVector3 v1, v2;
59 v1.SetMagThetaPhi(1, TMath::Pi()/2-fRawEl, fRawAz);
60 v2.SetMagThetaPhi(1, TMath::Pi()/2-fStarEl, fStarAz);
61
62 return v1.Angle(v2)*TMath::RadToDeg();
63 */
64
65 const Double_t del = fRawEl-fStarEl;
66 const Double_t daz = fRawAz-fStarAz;
67
68 /*
69 const Double_t dphi2 = daz/2.;
70 const Double_t cos2 = cos(dphi2)*cos(dphi2);
71 const Double_t sin2 = sin(dphi2)*sin(dphi2);
72 const Double_t d = cos(del)*cos2 - cos(fRawEl+fStarEl)*sin2;
73 */
74 const Double_t d = cos(del) - cos(fRawEl)*cos(fStarEl)*(1.-cos(daz));
75
76 if (err)
77 {
78 // Error of one pixel in the CCD
79 const Double_t e1 = 32./3600*TMath::DegToRad() * 0.5;
80
81 // Error of one SE unit
82 const Double_t e2 = 360./16384*TMath::DegToRad() * 0.5;
83
84 const Double_t e11 = sin(del)+cos(fRawEl)*sin(fStarEl)*(1-cos(daz));
85 const Double_t e12 = cos(fRawEl)*cos(fStarEl)*sin(daz);
86
87 const Double_t e21 = -sin(del)+sin(fRawEl)*cos(fStarEl)*(1-cos(daz));
88 const Double_t e22 = -cos(fRawEl)*cos(fStarEl)*sin(daz);
89
90 const Double_t err1 = sqrt(1-d*d);
91 const Double_t err2 = (e11*e11 + e12*e12)*e1*e1;
92 const Double_t err3 = (e21*e21 + e22*e22)*e2*e2;
93
94 *err = sqrt(err2+err3)/err1 * TMath::RadToDeg();
95 }
96
97 const Double_t dist = acos(d);
98 return dist * TMath::RadToDeg();
99 }
100
101istream &operator>>(istream &fin, TPointStar &set)
102{
103 TString str;
104 do
105 {
106 str.ReadLine(fin);
107 if (!fin)
108 return fin;
109 } while (str[0]=='#');
110
111 Float_t v[4], mag;
112 Int_t n = sscanf(str.Data(), "%f %f %f %f %*f %*f %*f %*f %*f %*f %f", v, v+1, v+2, v+3, &mag);
113 if (n<4)
114 {
115 gLog << err << "Read: ERROR - Not enough numbers" << endl;
116 return fin;
117 }
118 set.fMag = n<5 ? -25 : mag;
119
120 set.fStarAz = v[0]*TMath::DegToRad();
121 set.fStarEl = v[1]*TMath::DegToRad();
122
123 set.fRawAz = v[2]*TMath::DegToRad();
124 set.fRawEl = v[3]*TMath::DegToRad();
125
126 if (fin)
127 {
128 Double_t res, err;
129 res = set.GetResidual(&err);
130 gLog << inf << "Read: " << v[0] << " " << v[1] << " : " << v[2] << " " << v[3] << " : " << v[2]-v[0] << " " << v[3]-v[1] << " : " << res << " " << err << " " << err/res << endl;
131 }
132
133 return fin;
134}
135
136ostream &operator<<(ostream &out, TPointStar &set)
137{
138 out << Form("%8.3f", set.fStarAz*TMath::RadToDeg()) << " ";
139 out << Form("%8.3f", set.fStarEl*TMath::RadToDeg()) << " ";
140 out << Form("%8.3f", set.fRawAz*TMath::RadToDeg()) << " ";
141 out << Form("%8.1f", set.fRawEl*TMath::RadToDeg()) << " ";
142 out << set.fMag;
143
144 return out;
145}
Note: See TracBrowser for help on using the repository browser.