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

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