source: branches/Corsika7405Compatibility/macros/pointing.C

Last change on this file was 3957, checked in by tbretz, 21 years ago
*** empty log message ***
File size: 12.5 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Thomas Bretz, 5/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2004
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// pointing.C
28// ==========
29//
30// This macro is a demonstartion how check plots for a subsystem
31// (here the drive) can be made using Mars.
32//
33// In this case a merpped (root-) cc-report file is read in. The data
34// of the Drive branch is extracted using MReadReports and the
35// Stream-Id feature of MTaskList (second argument in AddToList).
36//
37// The output are plots showing (hopefully) the peformance of the system.
38//
39/////////////////////////////////////////////////////////////////////////////
40
41/*
42 class MGraph : public TGraph
43{
44public:
45 MGraph() : TGraph() {}
46 MGraph(Int_t n) : TGraph(n) {}
47 MGraph(Int_t n, const Int_t *x, const Int_t *y) : TGraph(n, x, y) {}
48 MGraph(Int_t n, const Float_t *x, const Float_t *y) : TGraph(n, x, y) {}
49 MGraph(Int_t n, const Double_t *x, const Double_t *y) : TGraph(n, x, y) {}
50 MGraph(const TVector &vx, const TVector &vy) : TGraph(vx, vy) {}
51 MGraph(const TVectorD &vx, const TVectorD &vy) : TGraph(vx, vy) {}
52 MGraph(const TH1 *h) : TGraph(h) {}
53 MGraph(const TF1 *f, Option_t *option="") : TGraph(f, option) {}
54 MGraph(const char *filename, const char *format="%lg %lg", Option_t *option="") : TGraph(filename, format, option) {}
55
56 void Paint(Option_t *o="")
57 {
58 if (!strstr(option,"POL") && !strstr(option,"pol"))
59 {
60 TGraph::Paint(o);
61 return;
62 }
63
64 //gPad->Range(-1.15, -1, 1.15, 1);
65
66 //gPad->Modified();
67 //gPad->Update();
68
69 TView *view = gPad->GetView();
70 if (!view)
71 {
72 cout << "No View!" << endl;
73 return;
74 }
75
76 TGraph gr;
77
78 Double_t *zd=g1->GetY();
79 Double_t *az=g2->GetY();
80
81 TMarker m;
82 m.SetMarkerStyle(kFullDotMedium);
83 m.SetMarkerColor(kRed);
84
85 for (int i=0; i<fN; i++)
86 {
87 const Double_t x = fX[i]/90;
88 const Double_t y = (fY[i]/180+1)*TMath::Pi();
89
90 Double_t r0[3] = { y*cos(x), y*sin(x), 0};
91 Double_t r1[3];
92
93 //gr.SetPoint(gr.GetN(), r0[0], r0[1]);
94
95 view->WCtoNDC(x, y);
96
97 m->PaintMarker(r1[0], r1[1]);
98 }
99 }
100};*/
101
102void pointing()
103{
104 MStatusDisplay *d = new MStatusDisplay;
105 d->SetLogStream(&gLog, kTRUE);
106
107 //
108 // Create a empty Parameter List and an empty Task List
109 // The tasklist is identified in the eventloop by its name
110 //
111 MParList plist;
112
113 MTaskList tlist;
114 plist.AddToList(&tlist);
115
116 //
117 // Now setup the tasks and tasklist:
118 // ---------------------------------
119 //
120
121 // First Task: Read file with image parameters
122 // (created with the star.C macro)
123 MReadReports read;
124 read.AddTree("Drive");
125 read.AddFile("data/2004_01_26_report.root");
126 read.AddToBranchList("MReportDrive.*");
127
128 MContinue tracking("MReportDrive.fState<3.5");
129
130 // Create a task which fills one histogram with the data
131 //MHVsTime h0("MReportDrive.fMjd");
132 MHVsTime h1("MReportDrive.fNominalZd");
133 MHVsTime h2("MReportDrive.fNominalAz");
134 MHVsTime h3("MReportDrive.fState");
135 MHVsTime h7("MReportDrive.GetAbsError*60");
136 //h0.SetName("Mjd");
137 h1.SetName("Zd");
138 h2.SetName("Az");
139 h3.SetName("State");
140 h7.SetName("ControlDeviation");
141
142 MH3 h4("MReportDrive.GetAbsError*60");
143 h4.SetName("DeltaH");
144
145 MH3 h5("MReportDrive.fNominalZd","MReportDrive.GetAbsError*60");
146 h5.SetName("DeltaHvsZd");
147
148 MBinning bins("BinningDeltaH");
149 bins.SetEdges(18, 0, 3.6);
150
151 MBinning bins2("BinningDeltaHvsZdX");
152 MBinning bins3("BinningDeltaHvsZdY");
153 bins2.SetEdges(90, 0, 90);
154 bins3.SetEdges(18, 0, 3.6);
155
156 plist.AddToList(&bins);
157 plist.AddToList(&bins2);
158 plist.AddToList(&bins3);
159
160 //MFillH fill0(&h0, "MTimeDrive");
161 MFillH fill1(&h1, "MTimeDrive");
162 MFillH fill2(&h2, "MTimeDrive");
163 MFillH fill3(&h3, "MTimeDrive");
164 MFillH fill7(&h7, "MTimeDrive");
165 MFillH fill4(&h4);
166 MFillH fill5(&h5);
167
168 //fill0.SetBit(MFillH::kDoNotDisplay);
169 fill1.SetBit(MFillH::kDoNotDisplay);
170 fill2.SetBit(MFillH::kDoNotDisplay);
171 fill3.SetBit(MFillH::kDoNotDisplay);
172 fill4.SetBit(MFillH::kDoNotDisplay);
173 fill5.SetBit(MFillH::kDoNotDisplay);
174 fill7.SetBit(MFillH::kDoNotDisplay);
175
176 // -----------------------------------------------
177
178 //
179 // Setup Task list
180 //
181 tlist.AddToList(&read);
182 tlist.AddToList(&fill3);
183 tlist.AddToList(&tracking);
184 //tlist.AddToList(&fill0);
185 tlist.AddToList(&fill1);
186 tlist.AddToList(&fill2);
187 tlist.AddToList(&fill4);
188 tlist.AddToList(&fill5);
189 tlist.AddToList(&fill7);
190
191 gStyle->SetOptStat(0);
192
193 //
194 // Create and setup the eventloop
195 //
196 MEvtLoop evtloop;
197 evtloop.SetParList(&plist);
198
199 //
200 // Execute your analysis
201 //
202 evtloop.SetDisplay(d);
203 if (!evtloop.Eventloop())
204 return;
205
206 tlist.PrintStatistics();
207
208
209 TGraph *g1 = h1.GetGraph();
210 TGraph *g2 = h2.GetGraph();
211
212 TCanvas &c = d->AddTab("Sky");
213 c.cd();
214
215 //Skyplot
216 TPad *p = new TPad("", "",0,0.32,0.5,1);
217 p->Draw();
218 p->cd();
219
220 gPad->SetTheta(-90);
221 gPad->SetPhi(90);
222 gPad->SetBorderMode(0);
223 gStyle->SetOptStat(0);
224
225 TH2F h("pol", "Telescope Tracking Positions on the Sky", 16, 0, 1, 9, 0, 1);
226 h.DrawClone("surf1pol");
227
228 gPad->Modified();
229 gPad->Update();
230
231 TView *view = gPad->GetView();
232 if (!view)
233 {
234 cout << "No View!" << endl;
235 return;
236 }
237
238 Double_t *zd=g1->GetY();
239 Double_t *az=g2->GetY();
240
241 Double_t old[2] = {0,0};
242
243 for (int i=0; i<g1->GetN(); i++)
244 {
245 az[i] += 180;
246 az[i] *= TMath::Pi()/180;
247
248 Double_t x[3] = { zd[i]*cos(az[i])/90, zd[i]*sin(az[i])/90, 0};
249 Double_t y[3];
250
251 view->WCtoNDC(x, y);
252
253 if (old[0]!=0 && old[1]!=1)
254 {
255 TLine *l = new TLine(y[0], y[1], old[0], old[1]);
256 l->SetLineColor(kBlue);
257 l->Draw();
258 }
259
260 TMarker *m = new TMarker(y[0], y[1], kFullDotMedium);
261 m->SetMarkerColor(i==g1->GetN()-1 ? kGreen : kRed);
262 m->Draw();
263
264 old[0] = y[0];
265 old[1] = y[1];
266 }
267
268 c.cd();
269
270
271/*
272 //MJD
273 p = new TPad("", "", 0.6, 0.66, 1, 1);
274 p->Draw();
275 p->cd();
276
277 MHVsTime *hvt=h0.DrawClone("nonew");
278 hvt->GetGraph()->SetMarkerStyle(kFullDotSmall);
279 hvt->GetGraph()->GetHistogram()->SetXTitle("Time");
280 hvt->GetGraph()->GetHistogram()->SetYTitle("");
281 hvt->GetGraph()->GetHistogram()->SetTitle("MJD vs. Time");
282 hvt->GetGraph()->GetHistogram()->SetStats(0);
283
284 c.cd();
285*/
286
287
288 //Histogram of Control Deviation
289 p = new TPad("", "", 0, 0, 0.5, 0.32);
290 p->Draw();
291 p->cd();
292
293 gPad->SetBorderMode(0);
294 //number of entries, mean, rms and number of overflows in the statusbox
295 gStyle->SetOptStat("emro");
296 //gStyle->SetStatFormat(".2g");
297
298 MH3 *mh3 = (MH3*)h4.DrawClone("nonew");
299
300 TAxis *axey = mh3->GetHist()->GetYaxis();
301 TAxis *axex = mh3->GetHist()->GetXaxis();
302 axey->SetLabelSize(0.05);
303 axex->SetLabelSize(0.05);
304 axex->SetTitleSize(0.05);
305 axex->SetTitleOffset(0.85);
306
307 mh3->GetHist()->SetXTitle("\\Delta [arcmin]");
308 mh3->GetHist()->SetYTitle("");
309 mh3->GetHist()->SetTitle("Control deviation of the motors");
310 mh3->GetHist()->SetStats(1);
311
312 //insert lines for 0.5, 1 and 2 SE
313 TLine ln;
314 ln.SetLineColor(kGreen);
315 ln.DrawLine(0.5*360*60/16384., 0, 0.5*360*60/16384., h4.GetHist()->GetMaximum());
316 ln.SetLineColor(kYellow);
317 ln.DrawLine(1.0*360*60/16384., 0, 1.0*360*60/16384., h4.GetHist()->GetMaximum());
318 ln.SetLineColor(kRed);
319 ln.DrawLine(2.0*360*60/16384., 0, 2.0*360*60/16384., h4.GetHist()->GetMaximum());
320
321 c.cd();
322
323
324 //Status of the Drive System vs Time
325 p = new TPad("", "", 0.5, 0.86, 1, 1);
326 p->Draw();
327 p->cd();
328 gPad->SetBorderMode(0);
329 p->SetGridx();
330
331 hvt = (MHVsTime*)h3.DrawClone("nonew");
332 hvt->GetGraph()->SetMarkerStyle(kFullDotSmall);
333
334 TH1 *hist = hvt->GetGraph()->GetHistogram();
335 TAxis *axey = hist->GetYaxis();
336 TAxis *axex = hist->GetXaxis();
337
338 hist->SetXTitle("Time");
339 hist->SetYTitle("");
340 hist->SetTitle("");//Drive Status vs. Time");
341 hist->SetStats(0);
342 hist->SetMinimum(-0.5);
343 hist->SetMaximum(4.5);
344 axey->Set(5, -0.5, 4.5);
345 axey->SetBinLabel(axey->FindFixBin(0), "Error");
346 axey->SetBinLabel(axey->FindFixBin(1), "Stopped");
347 axey->SetBinLabel(axey->FindFixBin(3), "Moving");
348 axey->SetBinLabel(axey->FindFixBin(4), "Tracking");
349 axey->SetLabelSize(0.15);
350 axex->SetLabelSize(0.08);
351 axex->SetTitleSize(0.09);
352 axex->SetTitleOffset(0.45);
353
354 c.cd();
355
356
357 //Zd vs Time
358 p = new TPad("", "", 0.5, 0.59, 1, 0.86);
359 p->Draw();
360 p->cd();
361 gPad->SetBorderMode(0);
362 p->SetGridx();
363
364 hvt = (MHVsTime*)h1.DrawClone("nonew");
365 hvt->GetGraph()->SetMarkerStyle(kFullDotSmall);
366
367 TH1 *hist = hvt->GetGraph()->GetHistogram();
368 if (hvt->GetGraph()->GetN())
369 {
370 hist->SetXTitle("Time");
371 hist->SetYTitle("Zd [\\circ]");
372 hist->SetTitle("Zd vs. Time");
373 hist->SetStats(0);
374 }
375
376 TAxis *axey = hist->GetYaxis();
377 TAxis *axex = hist->GetXaxis();
378 axey->SetLabelSize(0.05);
379 axey->SetTitleSize(0.06);
380 axey->SetTitleOffset(0.6);
381 axex->SetLabelSize(0.05);
382 axex->SetTitleSize(0.06);
383 axex->SetTitleOffset(0.85);
384
385 c.cd();
386
387
388 //Controll Deviation vs Time
389 p = new TPad("", "", 0.5, 0.32, 1, 0.59);
390 p->Draw();
391 p->cd();
392 gPad->SetBorderMode(0);
393 p->SetGridx();
394
395 hvt = (MHVsTime*)h7.DrawClone("nonew L");
396 TH1 *hist = hvt->GetGraph()->GetHistogram();
397 hist->SetAxisRange(-0.5, 10.5, "Y");
398 hist->SetXTitle("Time");
399 hist->SetYTitle("\\Delta [arcmin]");
400 hist->SetTitle("Control Deviation vs. Time");
401
402 TAxis *axey = hist->GetYaxis();
403 TAxis *axex = hist->GetXaxis();
404 axey->SetLabelSize(0.05);
405 axey->SetTitleSize(0.06);
406 axey->SetTitleOffset(0.5);
407 axex->SetLabelSize(0.05);
408 axex->SetTitleSize(0.06);
409 axex->SetTitleOffset(0.8);
410
411 //insert lines for 0.5, 1 and 2 SE
412 TLine ln;
413 ln.SetLineColor(kGreen);
414 ln.DrawLine(hist->GetXaxis()->GetXmin(), 0.5*360*60/16384., hist->GetXaxis()->GetXmax(), 0.5*360*60/16384.);
415 ln.SetLineColor(kYellow);
416 ln.DrawLine(hist->GetXaxis()->GetXmin(), 1.0*360*60/16384., hist->GetXaxis()->GetXmax(), 1.0*360*60/16384.);
417 ln.SetLineColor(kRed);
418 ln.DrawLine(hist->GetXaxis()->GetXmin(), 2.0*360*60/16384., hist->GetXaxis()->GetXmax(), 2.0*360*60/16384.);
419
420 c.cd();
421
422
423 //Controll Deviation vs Zd
424 p = new TPad("", "", 0.5, 0, 1, 0.32);
425 p->Draw();
426 p->cd();
427
428 gPad->SetBorderMode(0);
429 gStyle->SetOptStat(1110);
430 gStyle->SetStatFormat(".2g");
431
432 mh3 = (MH3*)h5.DrawClone("nonew");
433
434 TAxis *axey = mh3->GetHist()->GetYaxis();
435 TAxis *axex = mh3->GetHist()->GetXaxis();
436 axey->SetLabelSize(0.05);
437 axey->SetTitleSize(0.05);
438 axey->SetTitleOffset(0.7);
439 axex->SetLabelSize(0.05);
440 axex->SetTitleSize(0.05);
441 axex->SetTitleOffset(0.85);
442
443 mh3->GetHist()->SetYTitle("Zd [\\circ]");
444 mh3->GetHist()->SetXTitle("\\Delta [arcmin]");
445 mh3->GetHist()->SetTitle("Control deviation of the motors");
446 mh3->GetHist()->SetStats(1);
447 mh3->GetHist()->Draw("box");
448
449 //insert lines for 0.5, 1 and 2 SE
450 TLine ln;
451 ln.SetLineColor(kGreen);
452 ln.DrawLine(0, 0.5*360*60/16384., 90, 0.5*360*60/16384.);
453 ln.SetLineColor(kYellow);
454 ln.DrawLine(0, 1.0*360*60/16384., 90, 1.0*360*60/16384.);
455 ln.SetLineColor(kRed);
456 ln.DrawLine(0, 2.0*360*60/16384., 90, 2.0*360*60/16384.);
457
458
459 //d.SaveAsPS(2, "rep_files/CC_2003_11_27_22_31_59-new.ps");
460 //d.SaveAsRoot(2, "~/dev/Mars/rep_files/CC_2003_11_27_22_31_59-new.root");
461
462}
Note: See TracBrowser for help on using the repository browser.