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

Last change on this file since 8254 was 8113, checked in by tbretz, 19 years ago
*** empty log message ***
File size: 6.9 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 " << setw(23) << 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 (" << setw(3) << (Int_t)h->GetEntries() << "/";
112 cout << setw(3) << g[0].GetN() << ") " << flush;
113
114 return g[0].GetN();
115}
116
117struct Description_t
118{
119 const char *fName;
120 const char *fTitle;
121 const char *fFile;
122};
123
124const Int_t counts = 16;
125Description_t desc[counts] =
126{
127 // No good pointing model has been applied yet
128 // {"0404", "TPoints Residuals 4/2004" , "tpoint/tpoint0404.txt"},
129 // {"0405", "TPoints Residuals 5/2004" , "tpoint/tpoint0405.txt"},
130 // {"04081", "TPoints Residuals 8/2004-1" , "tpoint/tpoint0408-1.txt"},
131 {"04082", "TPoints Residuals 8/2004-2" , "tpoint/tpoint0408-2.txt"},
132 {"0409", "TPoints Residuals 9/2004" , "tpoint/tpoint0409.txt"},
133 // Culmination tests
134 // {"0410", "TPoints Residuals 10/2004" , "tpoint/tpoint0410.txt"},
135 {"0411", "TPoints Residuals 11/2004" , "tpoint/tpoint0411.txt"},
136 {"0412", "TPoints Residuals 12/2004" , "tpoint/tpoint0412.txt"},
137 // Worse pointing due to realignment of the mirror
138 {"0503", "TPoints Residuals 3/2005" , "tpoint/tpoint0503.txt"},
139 {"0504", "TPoints Residuals 4/2005" , "tpoint/tpoint0504.txt"},
140 {"05051", "TPoints Residuals 5/2005-1" , "tpoint/tpoint0505-1.txt"},
141 // Mirror alignment has been fixed
142 {"05052", "TPoints Residuals 5/2005-2" , "tpoint/tpoint0505-2.txt"},
143 // Fixes to pointing model due to fixing a screw
144 {"0506", "TPoints Residuals 6/2005" , "tpoint/tpoint0506.txt"},
145
146 {"0508", "TPoints Residuals 8/2005" , "tpoint/tpoint0508.txt"},
147 {"0509", "TPoints Residuals 9/2005" , "tpoint/tpoint0509.txt"},
148 // Quick-and-dirty mirror alignment
149 {"0510", "TPoints Residuals 10/2005" , "tpoint/tpoint0510.txt"},
150 // New mirror alignment after Tenerife meeting
151 {"0511", "TPoints Residuals 11/2005" , "tpoint/tpoint0511.txt"},
152
153 {"0512", "TPoints Residuals 12/2005" , "tpoint/tpoint0512.txt"},
154 {"0601", "TPoints Residuals 1/2006" , "tpoint/tpoint0601.txt"},
155 {"0603", "TPoints Residuals 3/2006" , "tpoint/tpoint0603.txt"}
156};
157
158void plot()
159{
160 TGraph g[3];
161
162 MBinning bins(10, 0, 0.2);
163
164 TArrayI n(counts);
165
166 TH1F hx[counts];
167 for (int i=0; i<counts; i++)
168 {
169 hx[i].SetNameTitle(desc[i].fName, desc[i].fTitle);
170 hx[i].SetDirectory(0);
171 bins.Apply(hx[i]);
172
173 cout << setw(2) << i << ": " << flush;
174 n[i] = fill(desc[i].fFile, g, &hx[i]);
175 cout << " Mean: " << setw(5) << setprecision(2) << hx[i].GetMean() << " deg" << hx[i].GetRMS()<<endl;
176 }
177
178 g[0].SetMarkerColor(kGreen);
179 g[1].SetMarkerColor(kMagenta);
180 g[2].SetMarkerColor(kBlack);
181 //g[2].SetLineColor(kBlack);
182 g[0].SetMarkerStyle(kFullDotMedium);
183 g[1].SetMarkerStyle(kFullDotMedium);
184 g[2].SetMarkerStyle(kFullDotLarge);
185
186 // --------- First Canvas ----------
187
188 new TCanvas("Vs Time", "");
189
190 TObject *obj[4];
191
192 for (int i=0; i<3; i++)
193 {
194 g[i].SetFillColor(kWhite);
195 g[i].SetLineColor(kWhite);
196 obj[i] = g[i].Clone();
197 obj[i]->SetBit(kCanDelete);
198 obj[i]->Draw(i==0?"AP":"P");
199 }
200
201 TLegend leg(0.905, 0.86, 0.99, 0.99);
202 leg.AddEntry(obj[0], " \\Delta\\theta");
203 leg.AddEntry(obj[1], " \\Delta\\phi");
204 leg.AddEntry(obj[2], " \\Delta");
205 leg.DrawClone()->SetBit(kCanDelete);
206
207 TLine l;
208 l.SetLineColor(kBlue);
209 for (int i=0; i<n.GetSize(); i++)
210 l.DrawLine(n[i], 0, n[i], 0.2);
211
212 // --------- Second Canvas ----------
213
214 new TCanvas("Distrib", "");
215
216 Double_t max=0;
217 for (int i=0; i<n.GetSize(); i++)
218 {
219 if (hx[i].GetEntries()==0)
220 {
221 cout << "Skip #" << i << endl;
222 continue;
223 }
224
225 hx[i].Scale(1./hx[i].GetEntries());
226 max = TMath::Max(max, hx[i].GetMaximum());
227 }
228 for (int i=0; i<n.GetSize(); i++)
229 {
230 hx[i].SetMaximum(max*1.05);
231 if (i<6)
232 hx[i].SetLineColor(kRed+i);
233 else
234 hx[i].SetMarkerStyle(kPlus+i-6);
235
236
237 }
238
239 for (int i=0; i<counts; i++)
240 hx[i].DrawCopy(i==0?"LP":"LPsame");
241
242 return;
243
244 for (int i=0; i<n.GetSize(); i++)
245 {
246 cout << "Mean: " << Form("%.3f +- %.3f", hx[i].GetMean(), hx[i].GetRMS());
247 cout << " (Overflows=" << hx[i].GetBinContent(hx[i].GetNbinsX()+1)*hx[i].GetEntries() << ")" << endl;
248 }
249}
Note: See TracBrowser for help on using the repository browser.