source: trunk/Cosy/tpoint/plot2.C@ 9966

Last change on this file since 9966 was 8984, checked in by tbretz, 17 years ago
*** empty log message ***
File size: 16.8 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 = 29+11+18+1;
104Description_t desc[counts] =
105{
106 /*
107 // 29. Apr. 2004 ~25800
108 // 5. Aug. 2004 ~32000
109 // 19. Aug. 2004 ~33530
110 // 7. Jun. 2005 ~57650
111 // 8. Jun. 2005
112 // 9. Jun. 2005 ~57860
113 // 12. Sep. 2005 ~68338
114 // 24. Nov. 2005 ~75562
115 // 17. Oct. 2006 ~103130
116 // 17. Jun. 2007 ~248193
117 */
118
119 // Culmination tests
120 //{"0411", "TPoints Residuals 11/2004" , "tpoint/tpoint0411.txt"},
121 //{"+0412", "TPoints Residuals 12/2004" , "tpoint/tpoint0412.txt"},
122
123 // 1: Worse pointing due to realignment of the mirror
124 {"0503", "TPoints Residuals 3/2005" , "tpoint/tpoint0503.txt"},
125
126 // New pointing model installed (29.4.2005)
127 // MIRROR MISALIGNMENT (WHEN?)
128 {"0504", "TPoints Residuals 4/2005" , "tpoint/tpoint0504.txt"},
129 {"+05051", "TPoints Residuals 5/2005-1" , "tpoint/tpoint0505-1.txt"},
130
131 // 2: Mirror alignment has been fixed
132 {"05052", "TPoints Residuals 5/2005-2" , "tpoint/tpoint0505-2.txt"},
133
134 // Pointing model changed due to fixing a screw
135 {"0506", "TPoints Residuals 6/2005" , "tpoint/tpoint0506.txt"},
136 // New pointing model applied (7.-9.6.2005)
137 {"0508", "TPoints Residuals 8/2005" , "tpoint/tpoint0508.txt"},
138
139 // New pointing model applied (12.9.2005)
140 {"0509", "TPoints Residuals 9/2005" , "tpoint/tpoint0509.txt"},
141 // Quick-and-dirty mirror alignment (only 4 TPoints)
142 {"+0510", "TPoints Residuals 10/2005" , "tpoint/tpoint0510.txt"},
143
144 // 3: New mirror alignment after Tenerife meeting
145 {"05111", "TPoints Residuals 11/2005-1" , "tpoint/tpoint0511-1.txt"},
146 //New pointing model installed (24.11.2005)
147 {"05112", "TPoints Residuals 11/2005-2" , "tpoint/tpoint0511-2.txt"},
148 {"+0512", "TPoints Residuals 12/2005" , "tpoint/tpoint0512.txt"},
149 {"+0601", "TPoints Residuals 1/2006" , "tpoint/tpoint0601.txt"},
150 {"+0603-1", "TPoints Residuals 3/2006-1" , "tpoint/tpoint0603-1.txt"},
151
152 // Changes to the mirror
153 {"0603-2", "TPoints Residuals 3/2006-2" , "tpoint/tpoint0603-2.txt"},
154 {"+0604", "TPoints Residuals 4/2006" , "tpoint/tpoint0604.txt"},
155 {"+0607", "TPoints Residuals 7/2006" , "tpoint/tpoint0607.txt"},
156 {"+0608", "TPoints Residuals 8/2006" , "tpoint/tpoint0608.txt"},
157 {"+0609", "TPoints Residuals 9/2006" , "tpoint/tpoint0609.txt"},
158 {"+0610", "TPoints Residuals 10/2006" , "tpoint/tpoint0610.txt"},
159
160 // 5: 6/10/17
161 {"0611", "TPoints Residuals 11/2006" , "tpoint/tpoint0611.txt"},
162 {"+0612", "TPoints Residuals 12/2006" , "tpoint/tpoint0612.txt"},
163 {"+0701", "TPoints Residuals 1/2007" , "tpoint/tpoint0701.txt"},
164 {"+0702", "TPoints Residuals 2/2007" , "tpoint/tpoint0702.txt"},
165 {"+0703", "TPoints Residuals 3/2007" , "tpoint/tpoint0703.txt"},
166
167 {"+0704", "TPoints Residuals 4/2007" , "tpoint/tpoint0704.txt"},
168 {"+0705", "TPoints Residuals 5/2007" , "tpoint/tpoint0705.txt"},
169 {"+0706", "TPoints Residuals 6/2007" , "tpoint/tpoint0706.txt"},
170 {"+07071", "TPoints Residuals 7/2007-1" , "tpoint/tpoint0707-1.txt"},
171 {"+07072", "TPoints Residuals 7/2007-2" , "tpoint/tpoint0707-2.txt"},
172 // 7/6/17
173
174 {"0708", "TPoints Residuals 7/8/23", "tpoint/2007_08_23_TPoints.txt"},
175 {"+0708", "TPoints Residuals 7/8/23", "tpoint/2007_08_24_TPoints.txt"},
176 {"+0708", "TPoints Residuals 7/8/23", "tpoint/2007_08_25_TPoints.txt"},
177 {"+0708", "TPoints Residuals 7/8/23", "tpoint/2007_08_26_TPoints.txt"},
178 {"+0708", "TPoints Residuals 7/8/23", "tpoint/2007_09_01_TPoints.txt"},
179 {"+0708", "TPoints Residuals 7/8/23", "tpoint/2007_09_03_TPoints.txt"},
180 {"+0708", "TPoints Residuals 7/8/23", "tpoint/2007_09_13_TPoints.txt"},
181 {"+0708", "TPoints Residuals 7/8/23", "tpoint/2007_09_15_TPoints.txt"},
182 {"+0708", "TPoints Residuals 7/8/23", "tpoint/2007_09_20_TPoints.txt"},
183 {"+0708", "TPoints Residuals 7/8/23", "tpoint/2007_10_04_TPoints.txt"},
184 {"+0708", "TPoints Residuals 7/8/23", "tpoint/2007_10_05_TPoints.txt"},
185
186 {"0801", "TPoints Residuals 7/8/23", "tpoint/2008_01_15_TPoints.txt"},
187 {"+0801", "TPoints Residuals 7/8/23", "tpoint/2008_01_18_TPoints.txt"},
188 {"+0801", "TPoints Residuals 7/8/23", "tpoint/2008_01_25_TPoints.txt"},
189 {"+0801", "TPoints Residuals 7/8/23", "tpoint/2008_01_26_TPoints.txt"},
190 {"+0801", "TPoints Residuals 7/8/23", "tpoint/2008_02_26_TPoints.txt"},
191 {"+0801", "TPoints Residuals 7/8/23", "tpoint/2008_03_16_TPoints.txt"},
192 {"+0801", "TPoints Residuals 7/8/23", "tpoint/2008_03_17_TPoints.txt"},
193 {"+0801", "TPoints Residuals 7/8/23", "tpoint/2008_03_18_TPoints.txt"},
194 {"+0801", "TPoints Residuals 7/8/23", "tpoint/2008_03_19_TPoints.txt"},
195 {"+0801", "TPoints Residuals 7/8/23", "tpoint/2008_03_25_TPoints.txt"},
196 {"+0801", "TPoints Residuals 7/8/23", "tpoint/2008_04_13_TPoints.txt"},
197 {"+0801", "TPoints Residuals 7/8/23", "tpoint/2008_04_14_TPoints.txt"},
198 {"+0801", "TPoints Residuals 7/8/23", "tpoint/2008_04_15_TPoints.txt"},
199 {"+0801", "TPoints Residuals 7/8/23", "tpoint/2008_04_16_TPoints.txt"},
200 {"+0801", "TPoints Residuals 7/8/23", "tpoint/2008_04_17_TPoints.txt"},
201 {"+0801", "TPoints Residuals 7/8/23", "tpoint/2008_04_23_TPoints.txt"},
202 {"+0801", "TPoints Residuals 7/8/23", "tpoint/2008_04_24_TPoints.txt"},
203 {"+0801", "TPoints Residuals 7/8/23", "tpoint/2008_04_25_TPoints.txt"},
204
205 // New pointing model 080611
206 {"0806", "TPoints Residuals 7/8/23", "tpoint/tpoint0806.txt"},
207 // New pointing model 080618
208};
209
210void plot2()
211{
212 TGraph g[3];
213
214 MBinning bins(100, 0, 0.2);
215
216 TArrayI n(counts);
217
218 TH1F hx[counts];
219 Int_t num = -1;
220 for (int i=0; i<counts; i++)
221 {
222 if (desc[i].fName[0]!='+')
223 {
224 num++;
225
226 hx[num].SetNameTitle(desc[i].fName, desc[i].fTitle);
227 hx[num].SetDirectory(0);
228 bins.Apply(hx[num]);
229 }
230
231 cout << setw(2) << num << ": " << flush;
232 n[num] = fill(desc[i].fFile, g, &hx[num]);
233 cout << " Mean: " << setw(5) << setprecision(2) << hx[num].GetMean() << " deg +/- " << hx[num].GetRMS()<< endl;
234 }
235
236 n.Set(++num);
237
238 g[0].SetMarkerColor(kGreen);
239 g[1].SetMarkerColor(kMagenta);
240 g[2].SetMarkerColor(kBlack);
241 //g[2].SetLineColor(kBlack);
242 g[0].SetMarkerStyle(kFullDotMedium);
243 g[1].SetMarkerStyle(kFullDotMedium);
244 g[2].SetMarkerStyle(kFullDotLarge);
245
246 // --------- First Canvas ----------
247
248 new TCanvas("Vs Time", "");
249
250 TObject *obj[4];
251
252 for (int i=0; i<3; i++)
253 {
254 g[i].SetFillColor(kWhite);
255 g[i].SetLineColor(kWhite);
256 obj[i] = g[i].Clone();
257 obj[i]->SetBit(kCanDelete);
258 obj[i]->Draw(i==0?"AP":"P");
259 }
260
261 TLegend leg(0.905, 0.86, 0.99, 0.99);
262 leg.AddEntry(obj[0], " \\Delta\\theta");
263 leg.AddEntry(obj[1], " \\Delta\\phi");
264 leg.AddEntry(obj[2], " \\Delta");
265 leg.DrawClone()->SetBit(kCanDelete);
266
267 TLine l;
268 l.SetLineColor(kBlue);
269 for (int i=0; i<n.GetSize(); i++)
270 l.DrawLine(n[i], 0, n[i], 0.2);
271
272 // --------- Second Canvas ----------
273
274 new TCanvas("Distrib", "");
275
276 Double_t max=0;
277 for (int i=0; i<n.GetSize(); i++)
278 {
279 if (hx[i].GetEntries()==0)
280 {
281 cout << "Skip #" << i << endl;
282 continue;
283 }
284
285 hx[i].Scale(1./hx[i].GetEntries());
286 max = TMath::Max(max, hx[i].GetMaximum());
287 }
288 for (int i=0; i<n.GetSize(); i++)
289 {
290 hx[i].SetMaximum(max*1.05);
291 if (i<6)
292 hx[i].SetLineColor(kRed+i);
293 else
294 hx[i].SetMarkerStyle(kPlus+i-6);
295
296
297 }
298
299 for (int i=0; i<counts; i++)
300 hx[i].DrawCopy(i==0?"LP":"LPsame");
301
302 Double_t time[] = {
303 MTime(2005, 3, 20).GetAxisTime(),
304 MTime(2005, 4, 29).GetAxisTime(),
305 MTime(2005, 5, 25).GetAxisTime(),
306 MTime(2005, 6, 8).GetAxisTime(),
307 MTime(2005, 8, 15).GetAxisTime(),
308 MTime(2005, 9, 12).GetAxisTime(),
309 MTime(2005, 11, 10).GetAxisTime(),
310 MTime(2005, 11, 24).GetAxisTime(),
311 MTime(2006, 3, 19).GetAxisTime(),
312 MTime(2006, 10, 17).GetAxisTime(),
313 MTime(2007, 7, 31).GetAxisTime(),
314 MTime(2008, 1, 14).GetAxisTime(),
315 MTime(2008, 6, 11).GetAxisTime(),
316 MTime(2008, 6, 18).GetAxisTime(),
317 };
318
319 TH1D histres[4];
320
321 MBinning bins;
322 bins.SetEdges(TArrayD(14, time));
323
324 bins.Apply(histres[0]);
325 bins.Apply(histres[1]);
326 bins.Apply(histres[2]);
327 bins.Apply(histres[3]);
328
329 TGraphAsymmErrors result[4];
330 TGraph resultm;
331 for (int i=0; i<n.GetSize(); i++)
332 {
333 cout << i+1 << " - Mean: " << Form("%.3f +- %.3f", hx[i].GetMean(), hx[i].GetRMS());
334 cout << " (Overflows=" << hx[i].GetBinContent(hx[i].GetNbinsX()+1)*hx[i].GetEntries() << ") " << (int)hx[i].GetEntries() << endl;
335
336 /*
337 TF1 fg("fg", "gaus", 0, 1);
338 hx[i].Fit(&fg);
339 result2.SetPoint(result.GetN(), result.GetN()+7, hx[i].GetMean());
340 result2.SetPointError(result.GetN()-1, 0, hx[i].GetRMS()/sqrt(hx[i].GetEntries()));
341 result.SetPoint(result.GetN(), result.GetN()+7, fg.GetParameter(1));
342 result.SetPointError(result.GetN()-1, 0, fg.GetParameter(2));
343 */
344
345 //Double_t q[4] = { MMath::GaussProb(0.5), MMath::GaussProb(1), MMath::GaussProb(2), MMath::GaussProb(3) };
346 Double_t q[4] = { 0.5, MMath::GaussProb(1), MMath::GaussProb(2), MMath::GaussProb(3) };
347 Double_t rc[4];
348 hx[i].GetQuantiles(4, rc, q);
349
350 for (int j=0; j<4; j++)
351 {
352 histres[j].SetBinContent(i+1, 60*rc[j]);
353
354 result[j].SetPoint(i, time[i], 60*rc[j]);
355 result[j].SetPointError(i, 0, 0, 60*rc[j], 0);
356
357 //result[j].SetPoint(result[j].GetN(), result[j].GetN()+1, 60*rc[j]);
358 //result[j].SetPointError(result[j].GetN()-1, 0, 0, 60*rc[j], 0);
359 }
360
361// result.SetPoint(result.GetN(), result.GetN()+1, rc[1]);
362// result2.SetPoint(result2.GetN(), result2.GetN()+1, rc[2]);
363
364 //result.SetPointError(result.GetN()-1, 0.5, 0);//rc[2]-rc[1]);
365 //result2.SetPointError(result.GetN()-1, 0.5, 0);//rc[2]-rc[1]);
366
367// result.SetPointError(result.GetN()-1, 0, 0, 0/*rc[1]-rc[0]*/, rc[2]-rc[1]);
368// result2.SetPointError(result.GetN()-1, 0, 0, 0/*rc[1]-rc[0]*/, rc[2]-rc[1]);
369
370// result2.SetPointError(result.GetN()-1, 0, 0, rc[2], 0);
371
372// resultm.SetPoint(resultm.GetN(), resultm.GetN()+1, 60*hx[i].GetMean());
373 resultm.SetPoint(resultm.GetN(), (time[i]+time[i+1])/2, 60*hx[i].GetMean());
374
375 }
376
377 new TCanvas;
378
379 gPad->SetBorderMode(0);
380 gPad->SetFrameBorderMode(0);
381 gPad->SetFillColor(kWhite);
382 gPad->SetRightMargin(0.01);
383 gPad->SetTopMargin(0.01);
384 gPad->SetLeftMargin(0.06);
385 gPad->SetGridy();
386
387 //Int_t col[] = { 12, 15, 17, 19 };
388 //Int_t col[] = { 12, 16, 18, 0 };
389 Int_t col[] = { 13, 16, 19, 0 };
390
391 TH1 *h = &histres[2];//result[3].GetHistogram();
392 h->SetXTitle("");
393 h->SetYTitle("Residual / arcmin");
394 h->SetBit(TH1::kNoStats);
395 h->GetXaxis()->CenterTitle();
396 h->GetYaxis()->CenterTitle();
397 h->GetYaxis()->SetTitleOffset(0.75);
398// h->GetXaxis()->SetTimeFormat("%m/%y %F1995-01-01 00:00:00 GMT");
399// h->GetXaxis()->SetTimeDisplay(1);
400
401 h->GetXaxis()->SetLabelColor(kWhite);
402
403 TLine line;
404
405 for (int j=2; j>=0; j--)
406 {
407 histres[j].SetMinimum(0);
408 histres[j].SetFillColor(col[j]);//12+2*j);
409
410 histres[j].DrawCopy(j==2?"":"same");
411
412 /*
413 //result[j].SetLineColor(kBlue);
414 //result[j].SetLineWidth(2);
415 //result[j].SetMarkerColor(kBlue);
416 result[j].SetMinimum(0);
417 //result2.SetMarkerStyle(kFullDotMedium);
418 //result[j].SetMarkerStyle(23);
419 result[j].SetFillColor(col[j]);//12+2*j);
420 result[j].DrawClone(j==3 ? "ABX" : "B"); // E3 B
421 */
422 }
423
424
425 resultm.SetMarkerStyle(20);
426 resultm.SetMarkerSize(0.8);
427 resultm.DrawClone("P");
428
429 line.DrawLine(time[0], 0, time[0], histres[2].GetMaximum());
430 for (int i=0; i<11; i++)
431 line.DrawLine(time[i+1], 0, time[i+1], histres[2].GetBinContent(i+1));
432
433 TText txt;
434 txt.SetTextSize(0.037);
435// txt.SetTextAngle(-45);
436 txt.SetTextAlign(23);
437 for (int m=4; m<13; m++)
438 {
439 Double_t monl = MTime(2005, m, 1,0).GetAxisTime();
440 Double_t mont = MTime(2005, m, 15,0).GetAxisTime();
441// txt.DrawText(mon, -0.12, Form("%02d/05", m));
442 txt.DrawText(mont, -0.10, Form("%d", m));
443 line.DrawLine(monl, -0.12, monl, 0.12);
444 }
445 for (int m=1; m<13; m++)
446 {
447 Double_t monl = MTime(2006, m, 1,0).GetAxisTime();
448 Double_t mont = MTime(2006, m, 15,0).GetAxisTime();
449// txt.DrawText(mon, -0.12, Form("%02d/06", m));
450 txt.DrawText(mont, -0.10, Form("%d", m));
451 line.DrawLine(monl, -0.12, monl, 0.12);
452 }
453 for (int m=1; m<13; m++)
454 {
455 Double_t monl = MTime(2007, m, 1,0).GetAxisTime();
456 Double_t mont = MTime(2007, m, 15,0).GetAxisTime();
457// txt.DrawText(mon, -0.12, Form("%02d/07", m));
458 txt.DrawText(mont, -0.10, Form("%d", m));
459 line.DrawLine(monl, -0.12, monl, 0.12);
460 }
461 for (int m=1; m<6; m++)
462 {
463 Double_t monl = MTime(2008, m, 1,0).GetAxisTime();
464 Double_t mont = MTime(2008, m, 15,0).GetAxisTime();
465// txt.DrawText(mon, -0.12, Form("%02d/07", m));
466 txt.DrawText(mont, -0.10, Form("%d", m));
467 line.DrawLine(monl, -0.12, monl, 0.12);
468 }
469
470 Double_t y6 = MTime(2006,1,1,0).GetAxisTime();
471 Double_t y7 = MTime(2007,1,1,0).GetAxisTime();
472 Double_t y8 = MTime(2008,1,1,0).GetAxisTime();
473 Double_t y9 = MTime(2008,5,1,0).GetAxisTime();
474
475 txt.SetTextSize(0.042);
476 txt.DrawText(y6-(y7-y6)/2, -0.6, "2005");
477 txt.DrawText((y6+y7)/2, -0.6, "2006");
478 txt.DrawText((y7+y8)/2, -0.6, "2007");
479 //txt.DrawText(y9+(y9-y8)/2, -0.6, "2008");
480
481 line.DrawLine(y6, -0.7, y6, 0.26);
482 line.DrawLine(y7, -0.7, y7, 0.26);
483 line.DrawLine(y8, -0.7, y8, 0.26);
484 line.SetLineStyle(3);
485 line.DrawLine(y6, 0.26, y6, 1.05*histres[2].GetMaximum());
486 line.DrawLine(y7, 0.26, y7, 1.05*histres[2].GetMaximum());
487 line.DrawLine(y8, 0.26, y8, 1.05*histres[2].GetMaximum());
488
489 line.SetLineColor(kBlue);
490 line.SetLineWidth(2);
491 line.SetLineStyle(kSolid);
492 line.DrawLine(time[0], 1*360/16384.*60, time[12], 1*360/16384.*60);
493 line.SetLineStyle(9);
494 line.DrawLine(time[0], 2*360/16384.*60, time[12], 2*360/16384.*60);
495 line.SetLineStyle(7);
496 line.DrawLine(time[0], 3*360/16384.*60, time[12], 3*360/16384.*60);
497
498 /* result.SetMinimum(-0.06);
499 result.SetMarkerStyle(kFullDotMedium);
500 result.DrawClone("LP");*/
501}
Note: See TracBrowser for help on using the repository browser.