source: trunk/MagicSoft/Cosy/plot.C@ 5245

Last change on this file since 5245 was 4725, checked in by tbretz, 21 years ago
*** empty log message ***
File size: 5.0 KB
Line 
1Double_t GetResidual(Double_t fRawEl, Double_t fRawAz,
2 Double_t fStarEl, Double_t fStarAz)
3{
4 fRawEl *= TMath::DegToRad();
5 fRawAz *= TMath::DegToRad();
6 fStarEl *= TMath::DegToRad();
7 fStarAz *= TMath::DegToRad();
8 Double_t del = fRawEl-fStarEl;
9 Double_t daz = fRawAz-fStarAz;
10 Double_t dphi2 = daz/2.;
11 Double_t cos2 = cos(dphi2)*cos(dphi2);
12 Double_t sin2 = sin(dphi2)*sin(dphi2);
13 Double_t d = cos(del)*cos2 - cos(fRawEl+fStarEl)*sin2;
14
15 Double_t dist = acos(d);
16
17 return dist * TMath::RadToDeg();
18}
19
20void DrawMarker(TVirtualPad *pad, Double_t r0, Double_t phi0, Double_t r1, Double_t phi1)
21 {
22 TView *view = pad->GetView();
23
24 if (!view)
25 {
26 cout << "No View!" << endl;
27 return;
28 }
29
30 TMarker mark0;
31 TMarker mark1;
32 mark0.SetMarkerStyle(kStar);
33 mark1.SetMarkerStyle(kStar);
34 mark1.SetMarkerColor(kRed);
35
36 r0 /= 90;
37 r1 /= 90;
38 phi0 *= TMath::DegToRad();
39 phi1 *= TMath::DegToRad();
40
41 Double_t x0[3] = { r0*cos(phi0), r0*sin(phi0), 0};
42 Double_t x1[3] = { r1*cos(phi1), r1*sin(phi1), 0};
43
44 mark0.DrawMarker(x0[0], x0[1]);
45 mark1.DrawMarker(x1[0], x1[1]);
46
47 return;
48 Double_t y0[3], y1[3];
49
50 view->WCtoNDC(x0, y0);
51 view->WCtoNDC(x1, y1);
52
53 mark0.DrawMarker(y0[0], y0[1]);
54 mark1.DrawMarker(y1[0], y1[1]);
55 }
56
57int fill(const char *fname, TGraph *g, TH1 *h)
58{
59/*
60 TH2F h2res1("Res2D1", " Dataset positions on the sky ", 32, 0, 360, 10, 0, 90);
61 h2res1.SetBit(TH1::kNoStats);
62 h2res1.DrawCopy("surf1pol");
63 gPad->Modified();
64 gPad->Update();
65 gPad->SetTheta(90);
66 gPad->SetPhi(-90);
67
68 DrawMarker(gPad, 45, 0, 0, 0);
69 gPad->Modified();
70 gPad->Update();
71
72 return;
73
74 */
75 ifstream fin(fname);
76
77 while (1)
78 {
79 TString str;
80 str.ReadLine(fin);
81 if (!fin)
82 break;
83
84 if (str.Contains("#"))
85 continue;
86
87 Float_t alt, az, dalt, daz, mjd;
88 sscanf(str.Data(), "%f %f %*f %*f %*f %*f %f %f %f",
89 &alt, &az, &dalt, &daz, &mjd);
90
91 if (dalt==0/* || GetResidual(alt, az, alt+dalt, az+daz)>0.1*/)
92 continue;
93
94 //cout << dalt << " " << daz << " " << GetResidual(alt, az, alt+dalt, az+daz) << endl;
95
96 mjd -= 53140.097505;
97
98 Double_t res = GetResidual(alt, az, alt-dalt, az-daz);
99
100 g[0].SetPoint(g[0].GetN(), g[0].GetN(), fabs(dalt));
101 g[1].SetPoint(g[1].GetN(), g[1].GetN(), fabs(daz));
102 g[2].SetPoint(g[2].GetN(), g[2].GetN(), res);
103
104 h->Fill(res);
105 }
106
107 return g[0].GetN();
108}
109
110void plot()
111{
112 TGraph g[3];
113
114 TH1F h0("0404", "TPoints Residuals 4/2004", 35, 0, 0.2);
115 TH1F h1("0405", "TPoints Residuals 5/2004", 35, 0, 0.2);
116 TH1F h2("0408", "TPoints Residuals 8/2004", 35, 0, 0.2);
117 TH1F h3("0409", "TPoints Residuals 9/2004", 35, 0, 0.2);
118
119 TH1 *h[4] = { &h0, &h1, &h2, &h3 };
120
121 TArrayI n(4);
122 n[0] = fill("tpoint0404.txt", g, &h0);
123 n[1] = fill("tpoint0405.txt", g, &h1);
124 n[2] = fill("tpoint0408.txt", g, &h2);
125 n[3] = fill("tpoint0409.txt", g, &h3);
126
127 for (int i=0; i<4; i++)
128 cout << "Overflows: " << Form("%4.0f", h[i]->GetBinContent(h[i]->GetNbinsX()+1)*h[i]->GetEntries()/h[i]->GetEntries()) << "/" << Form("%4.0f", h[i]->GetEntries()) << endl;
129
130 g[0].SetMarkerColor(kGreen);
131 g[1].SetMarkerColor(kMagenta);
132 g[2].SetMarkerColor(kBlack);
133 //g[2].SetLineColor(kBlack);
134 g[0].SetMarkerStyle(kFullDotMedium);
135 g[1].SetMarkerStyle(kFullDotMedium);
136 g[2].SetMarkerStyle(kFullDotLarge);
137
138 // --------- First Canvas ----------
139
140 new TCanvas("Vs Time", "");
141
142 TObject *obj[3];
143
144 for (int i=0; i<3; i++)
145 {
146 g[i].SetFillColor(kWhite);
147 g[i].SetLineColor(kWhite);
148 obj[i] = g[i].Clone();
149 obj[i]->SetBit(kCanDelete);
150 obj[i]->Draw(i==0?"AP":"P");
151 }
152
153 TLegend leg(0.905, 0.86, 0.99, 0.99);
154 leg.AddEntry(obj[0], " \\Delta\\theta");
155 leg.AddEntry(obj[1], " \\Delta\\phi");
156 leg.AddEntry(obj[2], " \\Delta");
157 leg.DrawClone()->SetBit(kCanDelete);
158
159 TLine l;
160 l.SetLineColor(kBlue);
161 for (int i=0; i<n.GetSize(); i++)
162 l.DrawLine(n[i], gPad->GetUymin(), n[i], gPad->GetUymax());
163
164 // --------- Second Canvas ----------
165
166 new TCanvas("Distrib", "");
167
168 Double_t max = 0;
169 for (int i=0; i<4; i++)
170 {
171 h[i]->Scale(1./h[i]->GetEntries());
172 max = TMath::Max(max, h[i]->GetMaximum());
173 }
174 for (int i=0; i<4; i++)
175 h[i]->SetMaximum(max*1.05);
176
177 h0.SetLineColor(kYellow);
178 h1.SetLineColor(kGreen);
179 h2.SetLineColor(kBlue);
180 //h3.SetLineColor(kYellow);
181
182 h3.DrawCopy();
183 h0.DrawCopy("same");
184 h1.DrawCopy("same");
185 h2.DrawCopy("same");
186 h3.DrawCopy("same");
187
188 for (int i=0; i<4; i++)
189 {
190 cout << "Mean: " << Form("%.3f +- %.3f", h[i]->GetMean(), h[i]->GetRMS());
191 cout << " (Overflows=" << h[i]->GetBinContent(h[i]->GetNbinsX()+1)*h[i]->GetEntries() << ")" << endl;
192 }
193}
Note: See TracBrowser for help on using the repository browser.