source: trunk/MagicSoft/Cosy/tpoint/plot.C@ 7217

Last change on this file since 7217 was 7147, checked in by tbretz, 19 years ago
*** empty log message ***
File size: 6.2 KB
Line 
1#include <iomanip>
2
3Double_t GetResidual(Double_t fRawEl, Double_t fRawAz,
4 Double_t fStarEl, Double_t fStarAz)
5{
6 fRawEl *= TMath::DegToRad();
7 fRawAz *= TMath::DegToRad();
8 fStarEl *= TMath::DegToRad();
9 fStarAz *= TMath::DegToRad();
10 Double_t del = fRawEl-fStarEl;
11 Double_t daz = fRawAz-fStarAz;
12 Double_t dphi2 = daz/2.;
13 Double_t cos2 = cos(dphi2)*cos(dphi2);
14 Double_t sin2 = sin(dphi2)*sin(dphi2);
15 Double_t d = cos(del)*cos2 - cos(fRawEl+fStarEl)*sin2;
16
17 Double_t dist = acos(d);
18
19 return dist * TMath::RadToDeg();
20}
21
22void DrawMarker(TVirtualPad *pad, Double_t r0, Double_t phi0, Double_t r1, Double_t phi1)
23 {
24 TView *view = pad->GetView();
25
26 if (!view)
27 {
28 cout << "No View!" << endl;
29 return;
30 }
31
32 TMarker mark0;
33 TMarker mark1;
34 mark0.SetMarkerStyle(kStar);
35 mark1.SetMarkerStyle(kStar);
36 mark1.SetMarkerColor(kRed);
37
38 r0 /= 90;
39 r1 /= 90;
40 phi0 *= TMath::DegToRad();
41 phi1 *= TMath::DegToRad();
42
43 Double_t x0[3] = { r0*cos(phi0), r0*sin(phi0), 0};
44 Double_t x1[3] = { r1*cos(phi1), r1*sin(phi1), 0};
45
46 mark0.DrawMarker(x0[0], x0[1]);
47 mark1.DrawMarker(x1[0], x1[1]);
48
49 return;
50 Double_t y0[3], y1[3];
51
52 view->WCtoNDC(x0, y0);
53 view->WCtoNDC(x1, y1);
54
55 mark0.DrawMarker(y0[0], y0[1]);
56 mark1.DrawMarker(y1[0], y1[1]);
57 }
58
59int fill(const char *fname, TGraph *g, TH1 *h)
60{
61/*
62 TH2F h2res1("Res2D1", " Dataset positions on the sky ", 32, 0, 360, 10, 0, 90);
63 h2res1.SetBit(TH1::kNoStats);
64 h2res1.DrawCopy("surf1pol");
65 gPad->Modified();
66 gPad->Update();
67 gPad->SetTheta(90);
68 gPad->SetPhi(-90);
69
70 DrawMarker(gPad, 45, 0, 0, 0);
71 gPad->Modified();
72 gPad->Update();
73
74 return;
75
76 */
77 ifstream fin(fname);
78
79 cout << "Reading from " << fname << "..." << flush;
80
81 while (1)
82 {
83 TString str;
84 str.ReadLine(fin);
85 if (!fin)
86 break;
87
88 if (str.Contains("#"))
89 continue;
90
91 Float_t alt, az, dalt, daz, mjd;
92 sscanf(str.Data(), "%f %f %*f %*f %*f %*f %f %f %f",
93 &alt, &az, &dalt, &daz, &mjd);
94
95 if (dalt==0/* || GetResidual(alt, az, alt+dalt, az+daz)>0.1*/)
96 continue;
97
98 //cout << dalt << " " << daz << " " << GetResidual(alt, az, alt+dalt, az+daz) << endl;
99
100 mjd -= 53140.097505;
101
102 Double_t res = GetResidual(alt, az, alt-dalt, az-daz);
103
104 g[0].SetPoint(g[0].GetN(), g[0].GetN(), fabs(dalt));
105 g[1].SetPoint(g[1].GetN(), g[1].GetN(), fabs(daz));
106 g[2].SetPoint(g[2].GetN(), g[2].GetN(), res);
107
108 h->Fill(res);
109 }
110
111 cout << "done (" << g[0].GetN() << ")" << endl;
112
113 return g[0].GetN();
114}
115
116struct Description_t
117{
118 const char *fName;
119 const char *fTitle;
120 const char *fFile;
121};
122
123const Int_t counts = 13-4;
124Description_t desc[counts] =
125{
126 // No good pointing model has been applied yet
127 // {"0404", "TPoints Residuals 4/2004" , "tpoint/tpoint0404.txt"},
128 // {"0405", "TPoints Residuals 5/2004" , "tpoint/tpoint0405.txt"},
129 // {"04081", "TPoints Residuals 8/2004-1" , "tpoint/tpoint0408-1.txt"},
130 {"04082", "TPoints Residuals 8/2004-2" , "tpoint/tpoint0408-2.txt"},
131 {"0409", "TPoints Residuals 9/2004" , "tpoint/tpoint0409.txt"},
132 // Culmination tests
133 // {"0410", "TPoints Residuals 10/2004" , "tpoint/tpoint0410.txt"},
134 {"0411", "TPoints Residuals 11/2004" , "tpoint/tpoint0411.txt"},
135 {"0412", "TPoints Residuals 12/2004" , "tpoint/tpoint0412.txt"},
136 // Worse pointing due to realignment of the mirror
137 {"0503", "TPoints Residuals 3/2005" , "tpoint/tpoint0503.txt"},
138 {"0504", "TPoints Residuals 4/2005" , "tpoint/tpoint0504.txt"},
139 {"05051", "TPoints Residuals 5/2004-1" , "tpoint/tpoint0505-1.txt"},
140 // Mirror alignment has been fixed
141 {"05052", "TPoints Residuals 5/2004-2" , "tpoint/tpoint0505-2.txt"},
142 // Fixes to pointing model due to fixing a screw
143 {"0506", "TPoints Residuals 6/2004" , "tpoint/tpoint0506.txt"}
144};
145
146void plot()
147{
148 TGraph g[3];
149
150 MBinning bins(10, 0, 0.2);
151
152 TArrayI n(counts);
153
154 TH1F hx[counts];
155 for (int i=0; i<counts; i++)
156 {
157 hx[i].SetNameTitle(desc[i].fName, desc[i].fTitle);
158 hx[i].SetDirectory(0);
159 bins.Apply(hx[i]);
160
161 cout << setw(2) << i << ": " << flush;
162 n[i] = fill(desc[i].fFile, g, &hx[i]);
163 }
164
165 g[0].SetMarkerColor(kGreen);
166 g[1].SetMarkerColor(kMagenta);
167 g[2].SetMarkerColor(kBlack);
168 //g[2].SetLineColor(kBlack);
169 g[0].SetMarkerStyle(kFullDotMedium);
170 g[1].SetMarkerStyle(kFullDotMedium);
171 g[2].SetMarkerStyle(kFullDotLarge);
172
173 // --------- First Canvas ----------
174
175 new TCanvas("Vs Time", "");
176
177 TObject *obj[4];
178
179 for (int i=0; i<3; i++)
180 {
181 g[i].SetFillColor(kWhite);
182 g[i].SetLineColor(kWhite);
183 obj[i] = g[i].Clone();
184 obj[i]->SetBit(kCanDelete);
185 obj[i]->Draw(i==0?"AP":"P");
186 }
187
188 TLegend leg(0.905, 0.86, 0.99, 0.99);
189 leg.AddEntry(obj[0], " \\Delta\\theta");
190 leg.AddEntry(obj[1], " \\Delta\\phi");
191 leg.AddEntry(obj[2], " \\Delta");
192 leg.DrawClone()->SetBit(kCanDelete);
193
194 TLine l;
195 l.SetLineColor(kBlue);
196 for (int i=0; i<n.GetSize(); i++)
197 l.DrawLine(n[i], 0, n[i], 0.2);
198
199 // --------- Second Canvas ----------
200
201 new TCanvas("Distrib", "");
202
203 Double_t max=0;
204 for (int i=0; i<n.GetSize(); i++)
205 {
206 if (hx[i].GetEntries()==0)
207 {
208 cout << "Skip #" << i << endl;
209 continue;
210 }
211
212 hx[i].Scale(1./hx[i].GetEntries());
213 max = TMath::Max(max, hx[i].GetMaximum());
214 }
215 for (int i=0; i<n.GetSize(); i++)
216 {
217 hx[i].SetMaximum(max*1.05);
218 if (i<6)
219 hx[i].SetLineColor(kRed+i);
220 else
221 hx[i].SetMarkerStyle(kPlus+i-6);
222
223
224 }
225
226 for (int i=0; i<counts; i++)
227 hx[i].DrawCopy(i==0?"LP":"LPsame");
228
229 return;
230
231 for (int i=0; i<n.GetSize(); i++)
232 {
233 cout << "Mean: " << Form("%.3f +- %.3f", hx[i].GetMean(), hx[i].GetRMS());
234 cout << " (Overflows=" << hx[i].GetBinContent(hx[i].GetNbinsX()+1)*hx[i].GetEntries() << ")" << endl;
235 }
236}
Note: See TracBrowser for help on using the repository browser.