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

Last change on this file since 9503 was 9496, checked in by tbretz, 16 years ago
*** empty log message ***
File size: 15.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 = 46;
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 {"+090512", "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_10/tpoint_20090609_232320.txt"},
138 {"+090512", "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_11/tpoint_20090611_011148.txt"},
139
140 // Something happened (reason unknown)
141
142 {"090611", "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_11/tpoint_20090611_023623.txt"},
143 {"+090611", "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_13/tpoint_20090613_033839.txt"},
144 {"+090611", "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_14/tpoint_20090614_021304.txt"},
145 {"+090611", "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_18/tpoint_20090618_041436.txt"},
146 {"+090611", "TPoints Residuals 11/2004" , "tpoint/m2/2009_07_12/tpoint_20090712_014259.txt"},
147 {"+090611", "TPoints Residuals 11/2004" , "tpoint/m2/2009_07_12/tpoint_20090712_024706.txt"},
148 {"+090611", "TPoints Residuals 11/2004" , "tpoint/m2/2009_07_13/tpoint_20090713_004244.txt"},
149 {"+090611", "TPoints Residuals 11/2004" , "tpoint/m2/2009_07_13/tpoint_20090713_025925.txt"},
150 {"+090611", "TPoints Residuals 11/2004" , "tpoint/m2/2009_07_14/tpoint_20090714_022322.txt"},
151 {"+090611", "TPoints Residuals 11/2004" , "tpoint/m2/2009_07_15/tpoint_20090715_021316.txt"},
152 {"+090611", "TPoints Residuals 08/2009" , "tpoint/m2/2009_08_02/tpoint_20090801_225305.txt"},
153 {"+090611", "TPoints Residuals 08/2009" , "tpoint/m2/2009_08_03/tpoint_20090802_224437.txt"},
154 {"+090611", "TPoints Residuals 08/2009" , "tpoint/m2/2009_08_04/tpoint_20090803_231853.txt"},
155 {"+090611", "TPoints Residuals 08/2009" , "tpoint/m2/2009_08_11/tpoint_20090810_234604.txt"},
156 {"+090611", "TPoints Residuals 08/2009" , "tpoint/m2/2009_08_11/tpoint_20090811_020844.txt"},
157 {"+090611", "TPoints Residuals 08/2009" , "tpoint/m2/2009_08_12/tpoint_20090812_001850.txt"},
158 {"+090611", "TPoints Residuals 08/2009" , "tpoint/m2/2009_08_13/tpoint_20090813_030805.txt"},
159
160 // ------ 09/08/17 new pointing model -------
161
162};
163
164void plot_m2()
165{
166 TGraph g[3];
167
168 MBinning bins(100, 0, 0.2);
169
170 TArrayI n(counts);
171
172 TH1F hx[counts];
173 Int_t num = -1;
174 for (int i=0; i<counts; i++)
175 {
176 if (desc[i].fName[0]!='+')
177 {
178 num++;
179
180 hx[num].SetNameTitle(desc[i].fName, desc[i].fTitle);
181 hx[num].SetDirectory(0);
182 bins.Apply(hx[num]);
183 }
184
185 cout << setw(2) << num << ": " << flush;
186 n[num] = fill(desc[i].fFile, g, &hx[num]);
187 cout << " Mean: " << setw(5) << setprecision(2) << hx[num].GetMean() << " deg +/- " << hx[num].GetRMS()<< endl;
188 }
189
190 n.Set(++num);
191
192 g[0].SetMarkerColor(kGreen);
193 g[1].SetMarkerColor(kMagenta);
194 g[2].SetMarkerColor(kBlack);
195 //g[2].SetLineColor(kBlack);
196 g[0].SetMarkerStyle(kFullDotMedium);
197 g[1].SetMarkerStyle(kFullDotMedium);
198 g[2].SetMarkerStyle(kFullDotLarge);
199
200 // --------- First Canvas ----------
201
202 new TCanvas("Vs Time", "");
203
204 TObject *obj[4];
205
206 for (int i=0; i<3; i++)
207 {
208 g[i].SetFillColor(kWhite);
209 g[i].SetLineColor(kWhite);
210 obj[i] = g[i].Clone();
211 obj[i]->SetBit(kCanDelete);
212 obj[i]->Draw(i==0?"AP":"P");
213 }
214
215 TLegend leg(0.905, 0.86, 0.99, 0.99);
216 leg.AddEntry(obj[0], " \\Delta\\theta");
217 leg.AddEntry(obj[1], " \\Delta\\phi");
218 leg.AddEntry(obj[2], " \\Delta");
219 leg.DrawClone()->SetBit(kCanDelete);
220
221 TLine l;
222 l.SetLineColor(kBlue);
223 for (int i=0; i<n.GetSize(); i++)
224 l.DrawLine(n[i], 0, n[i], 0.2);
225
226 // --------- Second Canvas ----------
227
228 new TCanvas("Distrib", "");
229
230 Double_t max=0;
231 for (int i=0; i<n.GetSize(); i++)
232 {
233 if (hx[i].GetEntries()==0)
234 {
235 cout << "Skip #" << i << endl;
236 continue;
237 }
238
239 hx[i].Scale(1./hx[i].GetEntries());
240 max = TMath::Max(max, hx[i].GetMaximum());
241 }
242 for (int i=0; i<n.GetSize(); i++)
243 {
244 hx[i].SetMaximum(max*1.05);
245 if (i<6)
246 hx[i].SetLineColor(kRed+i);
247 else
248 hx[i].SetMarkerStyle(kPlus+i-6);
249
250
251 }
252
253 for (int i=0; i<counts; i++)
254 hx[i].DrawCopy(i==0?"LP":"LPsame");
255
256 Double_t time[] = {
257 MTime(2009, 3, 30).GetAxisTime(),
258 MTime(2009, 4, 11).GetAxisTime(),
259 MTime(2009, 5, 11).GetAxisTime(),
260 MTime(2009, 6, 11).GetAxisTime(),
261 MTime(2009, 8, 17).GetAxisTime(),
262 MTime(2009, 12, 31).GetAxisTime(),
263 };
264
265 TH1D histres[4];
266
267 MBinning bins;
268 bins.SetEdges(TArrayD(6, time));
269
270 bins.Apply(histres[0]);
271 bins.Apply(histres[1]);
272 bins.Apply(histres[2]);
273 bins.Apply(histres[3]);
274
275 TGraphAsymmErrors result[4];
276 TGraph resultm;
277 for (int i=0; i<n.GetSize(); i++)
278 {
279 cout << i+1 << " - Mean: " << Form("%.4f +- %.4f", hx[i].GetMean(), hx[i].GetRMS());
280 cout << " (Overflows=" << hx[i].GetBinContent(hx[i].GetNbinsX()+1)*hx[i].GetEntries() << ") " << (int)hx[i].GetEntries() << endl;
281
282 /*
283 TF1 fg("fg", "gaus", 0, 1);
284 hx[i].Fit(&fg);
285 result2.SetPoint(result.GetN(), result.GetN()+7, hx[i].GetMean());
286 result2.SetPointError(result.GetN()-1, 0, hx[i].GetRMS()/sqrt(hx[i].GetEntries()));
287 result.SetPoint(result.GetN(), result.GetN()+7, fg.GetParameter(1));
288 result.SetPointError(result.GetN()-1, 0, fg.GetParameter(2));
289 */
290
291 //Double_t q[4] = { MMath::GaussProb(0.5), MMath::GaussProb(1), MMath::GaussProb(2), MMath::GaussProb(3) };
292 Double_t q[4] = { 0.5, MMath::GaussProb(1), MMath::GaussProb(2), MMath::GaussProb(3) };
293 Double_t rc[4];
294 hx[i].GetQuantiles(4, rc, q);
295
296 for (int j=0; j<4; j++)
297 {
298 histres[j].SetBinContent(i+1, 60*rc[j]);
299
300 result[j].SetPoint(i, time[i], 60*rc[j]);
301 result[j].SetPointError(i, 0, 0, 60*rc[j], 0);
302
303 //result[j].SetPoint(result[j].GetN(), result[j].GetN()+1, 60*rc[j]);
304 //result[j].SetPointError(result[j].GetN()-1, 0, 0, 60*rc[j], 0);
305 }
306
307// result.SetPoint(result.GetN(), result.GetN()+1, rc[1]);
308// result2.SetPoint(result2.GetN(), result2.GetN()+1, rc[2]);
309
310 //result.SetPointError(result.GetN()-1, 0.5, 0);//rc[2]-rc[1]);
311 //result2.SetPointError(result.GetN()-1, 0.5, 0);//rc[2]-rc[1]);
312
313// result.SetPointError(result.GetN()-1, 0, 0, 0/*rc[1]-rc[0]*/, rc[2]-rc[1]);
314// result2.SetPointError(result.GetN()-1, 0, 0, 0/*rc[1]-rc[0]*/, rc[2]-rc[1]);
315
316// result2.SetPointError(result.GetN()-1, 0, 0, rc[2], 0);
317
318// resultm.SetPoint(resultm.GetN(), resultm.GetN()+1, 60*hx[i].GetMean());
319 resultm.SetPoint(resultm.GetN(), (time[i]+time[i+1])/2, 60*hx[i].GetMean());
320
321 }
322
323 new TCanvas;
324
325 gPad->SetBorderMode(0);
326 gPad->SetFrameBorderMode(0);
327 gPad->SetFillColor(kWhite);
328 gPad->SetRightMargin(0.01);
329 gPad->SetTopMargin(0.01);
330 gPad->SetLeftMargin(0.06);
331 gPad->SetGridy();
332
333 //Int_t col[] = { 12, 15, 17, 19 };
334 //Int_t col[] = { 12, 16, 18, 0 };
335 Int_t col[] = { 13, 16, 19, 0 };
336
337 TH1 *h = &histres[2];//result[3].GetHistogram();
338 h->SetXTitle("");
339 h->SetYTitle("Residual / arcmin");
340 h->SetBit(TH1::kNoStats);
341 h->GetXaxis()->CenterTitle();
342 h->GetYaxis()->CenterTitle();
343 h->GetYaxis()->SetTitleOffset(0.75);
344// h->GetXaxis()->SetTimeFormat("%m/%y %F1995-01-01 00:00:00 GMT");
345// h->GetXaxis()->SetTimeDisplay(1);
346
347 h->GetXaxis()->SetLabelColor(kWhite);
348
349 TLine line;
350
351 for (int j=2; j>=0; j--)
352 {
353 histres[j].SetMinimum(0);
354 histres[j].SetFillColor(col[j]);//12+2*j);
355
356 histres[j].DrawCopy(j==2?"":"same");
357
358 /*
359 //result[j].SetLineColor(kBlue);
360 //result[j].SetLineWidth(2);
361 //result[j].SetMarkerColor(kBlue);
362 result[j].SetMinimum(0);
363 //result2.SetMarkerStyle(kFullDotMedium);
364 //result[j].SetMarkerStyle(23);
365 result[j].SetFillColor(col[j]);//12+2*j);
366 result[j].DrawClone(j==3 ? "ABX" : "B"); // E3 B
367 */
368 }
369
370
371 resultm.SetMarkerStyle(20);
372 resultm.SetMarkerSize(0.8);
373 resultm.DrawClone("P");
374
375 line.DrawLine(time[0], 0, time[0], histres[2].GetMaximum());
376 for (int i=0; i<bins.GetNumBins()-1; i++)
377 line.DrawLine(time[i+1], 0, time[i+1], histres[2].GetBinContent(i+1));
378
379 TText txt;
380 txt.SetTextSize(0.037);
381// txt.SetTextAngle(-45);
382 txt.SetTextAlign(23);
383 /*
384 for (int m=4; m<13; m++)
385 {
386 Double_t monl = MTime(2005, m, 1,0).GetAxisTime();
387 Double_t mont = MTime(2005, m, 15,0).GetAxisTime();
388// txt.DrawText(mon, -0.12, Form("%02d/05", m));
389 txt.DrawText(mont, -0.10, Form("%d", m));
390 line.DrawLine(monl, -0.12, monl, 0.12);
391 }
392 for (int m=1; m<13; m++)
393 {
394 Double_t monl = MTime(2006, m, 1,0).GetAxisTime();
395 Double_t mont = MTime(2006, m, 15,0).GetAxisTime();
396// txt.DrawText(mon, -0.12, Form("%02d/06", m));
397 txt.DrawText(mont, -0.10, Form("%d", m));
398 line.DrawLine(monl, -0.12, monl, 0.12);
399 }
400 for (int m=1; m<13; m++)
401 {
402 Double_t monl = MTime(2007, m, 1,0).GetAxisTime();
403 Double_t mont = MTime(2007, m, 15,0).GetAxisTime();
404// txt.DrawText(mon, -0.12, Form("%02d/07", m));
405 txt.DrawText(mont, -0.10, Form("%d", m));
406 line.DrawLine(monl, -0.12, monl, 0.12);
407 }*/
408 for (int m=4; m<13; m++)
409 {
410 Double_t monl = MTime(2009, m, 1,0).GetAxisTime();
411 Double_t mont = MTime(2009, m, 15,0).GetAxisTime();
412// txt.DrawText(mon, -0.12, Form("%02d/07", m));
413 txt.DrawText(mont, -0.10, Form("%d", m));
414 line.DrawLine(monl, -0.12, monl, 0.12);
415 }
416
417// Double_t y6 = MTime(2006,1,1,0).GetAxisTime();
418// Double_t y7 = MTime(2007,1,1,0).GetAxisTime();
419// Double_t y8 = MTime(2008,1,1,0).GetAxisTime();
420 Double_t y9 = MTime(2009,1,1,0).GetAxisTime();
421 Double_t y0 = MTime(2010,1,1,0).GetAxisTime();
422
423 txt.SetTextSize(0.042);
424// txt.DrawText(y6-(y7-y6)/2, -0.6, "2005");
425// txt.DrawText((y6+y7)/2, -0.6, "2006");
426// txt.DrawText((y7+y8)/2, -0.6, "2007");
427 txt.DrawText((y9+y0)/2, -0.6, "2009");
428
429// line.DrawLine(y6, -0.7, y6, 0.26);
430// line.DrawLine(y7, -0.7, y7, 0.26);
431// line.DrawLine(y8, -0.7, y8, 0.26);
432 line.DrawLine(y9, -0.7, y9, 0.26);
433 line.SetLineStyle(3);
434// line.DrawLine(y6, 0.26, y6, 1.05*histres[2].GetMaximum());
435// line.DrawLine(y7, 0.26, y7, 1.05*histres[2].GetMaximum());
436// line.DrawLine(y8, 0.26, y8, 1.05*histres[2].GetMaximum());
437 line.DrawLine(y9, 0.26, y9, 1.05*histres[2].GetMaximum());
438
439 line.SetLineColor(kBlue);
440 line.SetLineWidth(2);
441 line.SetLineStyle(kSolid);
442 line.DrawLine(time[0], 1*360/16384.*60, time[bins.GetNumBins()], 1*360/16384.*60);
443 line.SetLineStyle(9);
444 line.DrawLine(time[0], 2*360/16384.*60, time[bins.GetNumBins()], 2*360/16384.*60);
445 line.SetLineStyle(7);
446 line.DrawLine(time[0], 3*360/16384.*60, time[bins.GetNumBins()], 3*360/16384.*60);
447
448 /* result.SetMinimum(-0.06);
449 result.SetMarkerStyle(kFullDotMedium);
450 result.DrawClone("LP");*/
451}
Note: See TracBrowser for help on using the repository browser.