source: trunk/MagicSoft/Cosy/tpoint/plot_m2.C@ 9488

Last change on this file since 9488 was 9479, checked in by tbretz, 15 years ago
*** empty log message ***
File size: 13.9 KB
Line 
1#include <iomanip>
2
3void DrawMarker(TVirtualPad *pad, Double_t r0, Double_t phi0, Double_t r1, Double_t phi1)
4 {
5 TView *view = pad->GetView();
6
7 if (!view)
8 {
9 cout << "No View!" << endl;
10 return;
11 }
12
13 TMarker mark0;
14 TMarker mark1;
15 mark0.SetMarkerStyle(kStar);
16 mark1.SetMarkerStyle(kStar);
17 mark1.SetMarkerColor(kRed);
18
19 r0 /= 90;
20 r1 /= 90;
21 phi0 *= TMath::DegToRad();
22 phi1 *= TMath::DegToRad();
23
24 Double_t x0[3] = { r0*cos(phi0), r0*sin(phi0), 0};
25 Double_t x1[3] = { r1*cos(phi1), r1*sin(phi1), 0};
26
27 mark0.DrawMarker(x0[0], x0[1]);
28 mark1.DrawMarker(x1[0], x1[1]);
29
30 return;
31 Double_t y0[3], y1[3];
32
33 view->WCtoNDC(x0, y0);
34 view->WCtoNDC(x1, y1);
35
36 mark0.DrawMarker(y0[0], y0[1]);
37 mark1.DrawMarker(y1[0], y1[1]);
38 }
39
40int fill(const char *fname, TGraph *g, TH1 *h)
41{
42/*
43 TH2F h2res1("Res2D1", " Dataset positions on the sky ", 32, 0, 360, 10, 0, 90);
44 h2res1.SetBit(TH1::kNoStats);
45 h2res1.DrawCopy("surf1pol");
46 gPad->Modified();
47 gPad->Update();
48 gPad->SetTheta(90);
49 gPad->SetPhi(-90);
50
51 DrawMarker(gPad, 45, 0, 0, 0);
52 gPad->Modified();
53 gPad->Update();
54
55 return;
56
57 */
58 ifstream fin(fname);
59
60 cout << "Reading " << setw(23) << fname << "..." << flush;
61
62 while (1)
63 {
64 TString str;
65 str.ReadLine(fin);
66 if (!fin)
67 break;
68
69 if (str.Contains("#"))
70 continue;
71
72 Float_t alt, az, dalt, daz, mjd;
73 sscanf(str.Data(), "%f %f %*f %*f %*f %*f %f %f %f",
74 &az, &alt, &dalt, &daz, &mjd);
75
76 if (dalt==0/* || GetResidual(alt, az, alt+dalt, az+daz)>0.1*/)
77 continue;
78
79 mjd -= 53140.097505;
80
81 Double_t res = MAstro::GetDevAbs(90-(alt-dalt), -dalt, -daz);
82
83 g[0].SetPoint(g[0].GetN(), g[0].GetN(), fabs(dalt));
84 g[1].SetPoint(g[1].GetN(), g[1].GetN(), fabs(daz));
85 g[2].SetPoint(g[2].GetN(), g[2].GetN(), res);
86
87 h->Fill(res);
88 }
89
90 cout << "done (" << setw(3) << (Int_t)h->GetEntries() << "/";
91 cout << setw(4) << g[0].GetN() << ") " << flush;
92
93 return g[0].GetN();
94}
95
96struct Description_t
97{
98 const char *fName;
99 const char *fTitle;
100 const char *fFile;
101};
102
103const Int_t counts = 27;
104Description_t desc[counts] =
105{
106 {"090401", "TPoints Residuals 8/2004-2" , "tpoint/m2/first/tpoints_m2_1.txt"},
107 {"+090401", "TPoints Residuals 9/2004" , "tpoint/m2/first/tpoints_m2_2.txt"},
108 {"+090401", "TPoints Residuals 11/2004" , "tpoint/m2/first/tpoint_20090403_025841.txt"},
109 {"+090401", "TPoints Residuals 11/2004" , "tpoint/m2/first/tpoint_20090403_214619.txt"},
110 {"+090401", "TPoints Residuals 11/2004" , "tpoint/m2/first/tpoint_20090404_235955.txt"},
111 {"+090401", "TPoints Residuals 11/2004" , "tpoint/m2/first/tpoint_20090405_015920.txt"},
112 // New pointing model
113 {"090402", "TPoints Residuals 11/2004" , "tpoint/m2/first/tpoint_20090411_231731.txt"},
114 {"+090402", "TPoints Residuals 11/2004" , "tpoint/m2/first/tpoint_20090412_232919.txt"},
115 {"+090402", "TPoints Residuals 11/2004" , "tpoint/m2/first/tpoint_20090414_001746.txt"},
116 {"+090402", "TPoints Residuals 11/2004" , "tpoint/m2/first/tpoint_20090418_012742.txt"},
117 {"+090402", "TPoints Residuals 11/2004" , "tpoint/m2/first/tpoint_20090418_024412.txt"},
118 // New pointing model
119 {"090512", "TPoints Residuals 11/2004" , "tpoint/m2/2009_05_12/tpoint_20090512_011933.txt"},
120 {"+090512", "TPoints Residuals 11/2004" , "tpoint/m2/2009_05_12/tpoint_20090512_013112.txt"},
121 {"+090512", "TPoints Residuals 11/2004" , "tpoint/m2/2009_05_13/tpoint_20090512_210644.txt"},
122 {"+090512", "TPoints Residuals 11/2004" , "tpoint/m2/2009_05_13/tpoint_20090513_002758.txt"},
123 {"+090512", "TPoints Residuals 11/2004" , "tpoint/m2/2009_05_13/tpoint_20090513_025123.txt"},
124 {"+090512", "TPoints Residuals 11/2004" , "tpoint/m2/2009_05_13/tpoint_20090513_033504.txt"},
125 {"+090512", "TPoints Residuals 11/2004" , "tpoint/m2/2009_05_14/tpoint_20090514_013332.txt"},
126 // Just start of a new observation period
127 {"+090512", "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_01/tpoint_20090531_215139.txt"},
128 {"+090512", "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_01/tpoint_20090531_222548.txt"},
129 {"+090512", "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_02/tpoint_20090601_223001.txt"},
130 {"+090512", "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_03/tpoint_20090602_213508.txt"},
131 {"+090512", "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_03/tpoint_20090602_230944.txt"},
132 {"+090512", "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_03/tpoint_20090603_000616.txt"},
133 {"+090512", "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_03/tpoint_20090603_011935.txt"},
134 {"+090512", "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_04/tpoint_20090603_215841.txt"},
135 {"+090512", "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_05/tpoint_20090604_215941.txt"},
136};
137
138void plot_m2()
139{
140 TGraph g[3];
141
142 MBinning bins(100, 0, 0.2);
143
144 TArrayI n(counts);
145
146 TH1F hx[counts];
147 Int_t num = -1;
148 for (int i=0; i<counts; i++)
149 {
150 if (desc[i].fName[0]!='+')
151 {
152 num++;
153
154 hx[num].SetNameTitle(desc[i].fName, desc[i].fTitle);
155 hx[num].SetDirectory(0);
156 bins.Apply(hx[num]);
157 }
158
159 cout << setw(2) << num << ": " << flush;
160 n[num] = fill(desc[i].fFile, g, &hx[num]);
161 cout << " Mean: " << setw(5) << setprecision(2) << hx[num].GetMean() << " deg +/- " << hx[num].GetRMS()<< endl;
162 }
163
164 n.Set(++num);
165
166 g[0].SetMarkerColor(kGreen);
167 g[1].SetMarkerColor(kMagenta);
168 g[2].SetMarkerColor(kBlack);
169 //g[2].SetLineColor(kBlack);
170 g[0].SetMarkerStyle(kFullDotMedium);
171 g[1].SetMarkerStyle(kFullDotMedium);
172 g[2].SetMarkerStyle(kFullDotLarge);
173
174 // --------- First Canvas ----------
175
176 new TCanvas("Vs Time", "");
177
178 TObject *obj[4];
179
180 for (int i=0; i<3; i++)
181 {
182 g[i].SetFillColor(kWhite);
183 g[i].SetLineColor(kWhite);
184 obj[i] = g[i].Clone();
185 obj[i]->SetBit(kCanDelete);
186 obj[i]->Draw(i==0?"AP":"P");
187 }
188
189 TLegend leg(0.905, 0.86, 0.99, 0.99);
190 leg.AddEntry(obj[0], " \\Delta\\theta");
191 leg.AddEntry(obj[1], " \\Delta\\phi");
192 leg.AddEntry(obj[2], " \\Delta");
193 leg.DrawClone()->SetBit(kCanDelete);
194
195 TLine l;
196 l.SetLineColor(kBlue);
197 for (int i=0; i<n.GetSize(); i++)
198 l.DrawLine(n[i], 0, n[i], 0.2);
199
200 // --------- Second Canvas ----------
201
202 new TCanvas("Distrib", "");
203
204 Double_t max=0;
205 for (int i=0; i<n.GetSize(); i++)
206 {
207 if (hx[i].GetEntries()==0)
208 {
209 cout << "Skip #" << i << endl;
210 continue;
211 }
212
213 hx[i].Scale(1./hx[i].GetEntries());
214 max = TMath::Max(max, hx[i].GetMaximum());
215 }
216 for (int i=0; i<n.GetSize(); i++)
217 {
218 hx[i].SetMaximum(max*1.05);
219 if (i<6)
220 hx[i].SetLineColor(kRed+i);
221 else
222 hx[i].SetMarkerStyle(kPlus+i-6);
223
224
225 }
226
227 for (int i=0; i<counts; i++)
228 hx[i].DrawCopy(i==0?"LP":"LPsame");
229
230 Double_t time[] = {
231 MTime(2009, 3, 30).GetAxisTime(),
232 MTime(2009, 4, 11).GetAxisTime(),
233 MTime(2009, 5, 11).GetAxisTime(),
234 MTime(2009, 5, 31).GetAxisTime(),
235 MTime(2009, 6, 6).GetAxisTime(),
236 };
237
238 TH1D histres[4];
239
240 MBinning bins;
241 bins.SetEdges(TArrayD(5, time));
242
243 bins.Apply(histres[0]);
244 bins.Apply(histres[1]);
245 bins.Apply(histres[2]);
246 bins.Apply(histres[3]);
247
248 TGraphAsymmErrors result[4];
249 TGraph resultm;
250 for (int i=0; i<n.GetSize(); i++)
251 {
252 cout << i+1 << " - Mean: " << Form("%.4f +- %.4f", hx[i].GetMean(), hx[i].GetRMS());
253 cout << " (Overflows=" << hx[i].GetBinContent(hx[i].GetNbinsX()+1)*hx[i].GetEntries() << ") " << (int)hx[i].GetEntries() << endl;
254
255 /*
256 TF1 fg("fg", "gaus", 0, 1);
257 hx[i].Fit(&fg);
258 result2.SetPoint(result.GetN(), result.GetN()+7, hx[i].GetMean());
259 result2.SetPointError(result.GetN()-1, 0, hx[i].GetRMS()/sqrt(hx[i].GetEntries()));
260 result.SetPoint(result.GetN(), result.GetN()+7, fg.GetParameter(1));
261 result.SetPointError(result.GetN()-1, 0, fg.GetParameter(2));
262 */
263
264 //Double_t q[4] = { MMath::GaussProb(0.5), MMath::GaussProb(1), MMath::GaussProb(2), MMath::GaussProb(3) };
265 Double_t q[4] = { 0.5, MMath::GaussProb(1), MMath::GaussProb(2), MMath::GaussProb(3) };
266 Double_t rc[4];
267 hx[i].GetQuantiles(4, rc, q);
268
269 for (int j=0; j<4; j++)
270 {
271 histres[j].SetBinContent(i+1, 60*rc[j]);
272
273 result[j].SetPoint(i, time[i], 60*rc[j]);
274 result[j].SetPointError(i, 0, 0, 60*rc[j], 0);
275
276 //result[j].SetPoint(result[j].GetN(), result[j].GetN()+1, 60*rc[j]);
277 //result[j].SetPointError(result[j].GetN()-1, 0, 0, 60*rc[j], 0);
278 }
279
280// result.SetPoint(result.GetN(), result.GetN()+1, rc[1]);
281// result2.SetPoint(result2.GetN(), result2.GetN()+1, rc[2]);
282
283 //result.SetPointError(result.GetN()-1, 0.5, 0);//rc[2]-rc[1]);
284 //result2.SetPointError(result.GetN()-1, 0.5, 0);//rc[2]-rc[1]);
285
286// result.SetPointError(result.GetN()-1, 0, 0, 0/*rc[1]-rc[0]*/, rc[2]-rc[1]);
287// result2.SetPointError(result.GetN()-1, 0, 0, 0/*rc[1]-rc[0]*/, rc[2]-rc[1]);
288
289// result2.SetPointError(result.GetN()-1, 0, 0, rc[2], 0);
290
291// resultm.SetPoint(resultm.GetN(), resultm.GetN()+1, 60*hx[i].GetMean());
292 resultm.SetPoint(resultm.GetN(), (time[i]+time[i+1])/2, 60*hx[i].GetMean());
293
294 }
295
296 new TCanvas;
297
298 gPad->SetBorderMode(0);
299 gPad->SetFrameBorderMode(0);
300 gPad->SetFillColor(kWhite);
301 gPad->SetRightMargin(0.01);
302 gPad->SetTopMargin(0.01);
303 gPad->SetLeftMargin(0.06);
304 gPad->SetGridy();
305
306 //Int_t col[] = { 12, 15, 17, 19 };
307 //Int_t col[] = { 12, 16, 18, 0 };
308 Int_t col[] = { 13, 16, 19, 0 };
309
310 TH1 *h = &histres[2];//result[3].GetHistogram();
311 h->SetXTitle("");
312 h->SetYTitle("Residual / arcmin");
313 h->SetBit(TH1::kNoStats);
314 h->GetXaxis()->CenterTitle();
315 h->GetYaxis()->CenterTitle();
316 h->GetYaxis()->SetTitleOffset(0.75);
317// h->GetXaxis()->SetTimeFormat("%m/%y %F1995-01-01 00:00:00 GMT");
318// h->GetXaxis()->SetTimeDisplay(1);
319
320 h->GetXaxis()->SetLabelColor(kWhite);
321
322 TLine line;
323
324 for (int j=2; j>=0; j--)
325 {
326 histres[j].SetMinimum(0);
327 histres[j].SetFillColor(col[j]);//12+2*j);
328
329 histres[j].DrawCopy(j==2?"":"same");
330
331 /*
332 //result[j].SetLineColor(kBlue);
333 //result[j].SetLineWidth(2);
334 //result[j].SetMarkerColor(kBlue);
335 result[j].SetMinimum(0);
336 //result2.SetMarkerStyle(kFullDotMedium);
337 //result[j].SetMarkerStyle(23);
338 result[j].SetFillColor(col[j]);//12+2*j);
339 result[j].DrawClone(j==3 ? "ABX" : "B"); // E3 B
340 */
341 }
342
343
344 resultm.SetMarkerStyle(20);
345 resultm.SetMarkerSize(0.8);
346 resultm.DrawClone("P");
347
348 line.DrawLine(time[0], 0, time[0], histres[2].GetMaximum());
349 for (int i=0; i<bins.GetNumBins(); i++)
350 line.DrawLine(time[i+1], 0, time[i+1], histres[2].GetBinContent(i+1));
351
352 TText txt;
353 txt.SetTextSize(0.037);
354// txt.SetTextAngle(-45);
355 txt.SetTextAlign(23);
356 /*
357 for (int m=4; m<13; m++)
358 {
359 Double_t monl = MTime(2005, m, 1,0).GetAxisTime();
360 Double_t mont = MTime(2005, m, 15,0).GetAxisTime();
361// txt.DrawText(mon, -0.12, Form("%02d/05", m));
362 txt.DrawText(mont, -0.10, Form("%d", m));
363 line.DrawLine(monl, -0.12, monl, 0.12);
364 }
365 for (int m=1; m<13; m++)
366 {
367 Double_t monl = MTime(2006, m, 1,0).GetAxisTime();
368 Double_t mont = MTime(2006, m, 15,0).GetAxisTime();
369// txt.DrawText(mon, -0.12, Form("%02d/06", m));
370 txt.DrawText(mont, -0.10, Form("%d", m));
371 line.DrawLine(monl, -0.12, monl, 0.12);
372 }
373 for (int m=1; m<13; m++)
374 {
375 Double_t monl = MTime(2007, m, 1,0).GetAxisTime();
376 Double_t mont = MTime(2007, m, 15,0).GetAxisTime();
377// txt.DrawText(mon, -0.12, Form("%02d/07", m));
378 txt.DrawText(mont, -0.10, Form("%d", m));
379 line.DrawLine(monl, -0.12, monl, 0.12);
380 }*/
381 for (int m=1; m<13; m++)
382 {
383 Double_t monl = MTime(2009, m, 1,0).GetAxisTime();
384 Double_t mont = MTime(2009, m, 15,0).GetAxisTime();
385// txt.DrawText(mon, -0.12, Form("%02d/07", m));
386 txt.DrawText(mont, -0.10, Form("%d", m));
387 line.DrawLine(monl, -0.12, monl, 0.12);
388 }
389
390// Double_t y6 = MTime(2006,1,1,0).GetAxisTime();
391// Double_t y7 = MTime(2007,1,1,0).GetAxisTime();
392// Double_t y8 = MTime(2008,1,1,0).GetAxisTime();
393 Double_t y9 = MTime(2009,1,1,0).GetAxisTime();
394 Double_t y0 = MTime(2010,1,1,0).GetAxisTime();
395
396 txt.SetTextSize(0.042);
397// txt.DrawText(y6-(y7-y6)/2, -0.6, "2005");
398// txt.DrawText((y6+y7)/2, -0.6, "2006");
399// txt.DrawText((y7+y8)/2, -0.6, "2007");
400 txt.DrawText((y9+y0)/2, -0.6, "2009");
401
402// line.DrawLine(y6, -0.7, y6, 0.26);
403// line.DrawLine(y7, -0.7, y7, 0.26);
404// line.DrawLine(y8, -0.7, y8, 0.26);
405 line.DrawLine(y9, -0.7, y9, 0.26);
406 line.SetLineStyle(3);
407// line.DrawLine(y6, 0.26, y6, 1.05*histres[2].GetMaximum());
408// line.DrawLine(y7, 0.26, y7, 1.05*histres[2].GetMaximum());
409// line.DrawLine(y8, 0.26, y8, 1.05*histres[2].GetMaximum());
410 line.DrawLine(y9, 0.26, y9, 1.05*histres[2].GetMaximum());
411
412 line.SetLineColor(kBlue);
413 line.SetLineWidth(2);
414 line.SetLineStyle(kSolid);
415 line.DrawLine(time[0], 1*360/16384.*60, time[bins.GetNumBins()], 1*360/16384.*60);
416 line.SetLineStyle(9);
417 line.DrawLine(time[0], 2*360/16384.*60, time[bins.GetNumBins()], 2*360/16384.*60);
418 line.SetLineStyle(7);
419 line.DrawLine(time[0], 3*360/16384.*60, time[bins.GetNumBins()], 3*360/16384.*60);
420
421 /* result.SetMinimum(-0.06);
422 result.SetMarkerStyle(kFullDotMedium);
423 result.DrawClone("LP");*/
424}
Note: See TracBrowser for help on using the repository browser.