source: trunk/MagicSoft/Cosy/main/MCaos.cc@ 4051

Last change on this file since 4051 was 2615, checked in by tbretz, 22 years ago
*** empty log message ***
File size: 9.8 KB
Line 
1#include "MCaos.h"
2
3#include <iostream>
4#include <iomanip>
5
6#include <TSystem.h>
7#include <TFile.h>
8#include <TTree.h>
9#include <TH1.h>
10#include <TH2.h>
11#include <TGraph.h>
12#include <TCanvas.h>
13
14#include "MTime.h"
15
16#include "Led.h"
17#include "Ring.h"
18#include "Rings.h"
19#include "FilterLed.h"
20
21#include "coord.h"
22
23void MCaos::ReadResources(const char *name="leds.txt")
24{
25 ifstream fin(name);
26 if (!fin)
27 {
28 cout << "ERROR - Cannot open " << name << endl;
29 return;
30 }
31
32 fPositions.Clear();
33
34 cout << " Reading " << name << ":" << endl;
35 cout << "------------------------------" << endl;
36
37 while (1)
38 {
39 Double_t px, py, ox, oy;
40 fin >> px >> py >> ox >> oy;
41 if (!fin)
42 break;
43
44 cout << " Led #" << fPositions.GetEntriesFast() << ": ";
45 cout << setw(3) << px << " ";
46 cout << setw(3) << py << " (";
47 cout << setw(3) << ox << ", ";
48 cout << setw(3) << oy << ")" << endl;
49 AddPosition(px, py, ox, oy);
50 }
51 cout << "Found " << fPositions.GetEntriesFast() << " leds." << endl;
52}
53
54void MCaos::OpenFile()
55{
56 int i=0;
57 char name[100];
58 while (1)
59 {
60 sprintf(name, "data/data%03d.root", i++);
61 if (gSystem->AccessPathName(name, kFileExists))
62 break;
63 }
64
65 if (fFile)
66 delete fFile;
67
68 fFile = new TFile(name, "RECREATE");
69
70 if (!fFile->IsOpen())
71 {
72 delete fFile;
73 fFile = NULL;
74
75 cout << "Error: Cannot open file '" << name << "'" << endl;
76 }
77
78 TTree *tree = new TTree("Data", "Real CaOs Data");
79
80 fLeds = new Leds;
81 fEvtTime = 0;
82
83 tree->Branch("Leds", "TClonesArray", &fLeds);
84 tree->Branch("ZenithDist.", &fZenithDist, "fZenithDist/D");
85 tree->Branch("Azimuth.", &fAzimuth, "fAzimuth/D");
86 tree->Branch("EvtTime.", &fEvtTime, "fEvtTime/D");
87
88 cout << "Root file '" << name << "' open." << endl;
89}
90
91void MCaos::CloseFile()
92{
93 if (!fFile)
94 return;
95
96 const TString name = fFile->GetName();
97 const Double_t n = ((TTree*)fFile->Get("Data"))->GetEntries();
98
99 fFile->Write();
100 delete fFile;
101 fFile = NULL;
102
103 cout << "Root file closed (n=" << n << ")" << endl;
104
105 if (n<1)
106 {
107 gSystem->Unlink(name);
108 cout << "Root file deleted - no entries." << endl;
109 }
110}
111
112void MCaos::InitHistograms()
113{
114 if (fHistpr)
115 return;
116
117 Rings r;
118 r.CalcRings(fPositions);
119
120 const Ring &c = r.GetCenter();
121
122 Double_t xmin = c.GetX()-50;
123 Double_t xmax = c.GetX()+50;
124
125 Double_t ymin = c.GetY()-50;
126 Double_t ymax = c.GetY()+50;
127
128 Double_t rmin = c.GetR()-50;
129 Double_t rmax = c.GetR()+50;
130
131 Int_t xbin = 1001;
132 Int_t ybin = 1001;
133 Int_t rbin = 1001;
134
135 const Int_t sz = 50;
136 fHistled = new TH2F*[fPositions.GetEntriesFast()];
137 fHistw = new TH1F*[fPositions.GetEntriesFast()];
138 for (int i=0; i<fPositions.GetEntriesFast(); i++)
139 {
140 TString name = "LED";
141 TString title = "LED #";
142
143 name += i;
144 title += i;
145
146 const Led &p = fPositions(i);
147 fHistled[i] = new TH2F(name, title,
148 20*sz+1, p.GetX()-sz, p.GetX()+sz,
149 20*sz+1, p.GetY()-sz, p.GetY()+sz);
150 fHistled[i]->SetXTitle("x [pix]");
151 fHistled[i]->SetYTitle("counts");
152
153 name = "Angle";
154 title = "Angle of the Led #";
155
156 name += i;
157 title += i;
158
159 fHistw[i] = new TH1F;
160 fHistw[i]->SetNameTitle(name, title);
161 fHistw[i]->SetBins(721, -180.5, 180.5);
162 fHistw[i]->SetXTitle("\\Phi [deg]");
163 fHistw[i]->SetYTitle("counts");
164 }
165
166 fHistallw = new TH1F;
167 fHistallw->SetNameTitle("allw","Rotation angel");
168 fHistallw->SetBins(26, -25, 25);
169 fHistallw->SetXTitle("\\Phi [arcmin]");
170 fHistallw->SetYTitle("counts");
171
172 fHistprxpry = new TH2F;
173 fHistprxpry->SetNameTitle("prx und pry","x- and y-coordniate of the ring-center");
174 fHistprxpry->SetBins(xbin, xmin, xmax, ybin, ymin, ymax);
175 fHistprxpry->SetXTitle("x [mm]");
176 fHistprxpry->SetYTitle("y [mm]");
177 fHistprxpry->SetZTitle("counts");
178
179 fGraphprx = new TGraph;
180 fGraphprx->SetTitle("time-developement of the x-coordinate of the ring-center");
181
182 fGraphpry = new TGraph;
183 fGraphpry->SetTitle("time-developement of the y-coordinate of the ring-center");
184
185 fGraphw = new TGraph*[fPositions.GetEntriesFast()];
186 for (int i=0; i<fPositions.GetEntriesFast(); i++)
187 {
188 TString title = "Time-developement of the angle ";
189 title += i;
190
191 fGraphw[i] = new TGraph;
192 fGraphw[i]->SetTitle(title);
193 }
194
195 fHistpr = new TH1F("pr","Radius of the ring", rbin, rmin, rmax);
196 fHistpr->SetXTitle("r [pix]");
197 fHistpr->SetYTitle("counts");
198}
199
200void MCaos::DeleteHistograms()
201{
202 TH1F *dummy = fHistpr;
203 fHistpr=NULL;
204
205 if (!dummy)
206 return;
207
208 delete dummy;
209 delete fHistprxpry;
210 delete fHistallw;
211 delete fGraphprx;
212 delete fGraphpry;
213
214 for (int i=0; i<6; i++)
215 {
216 delete fHistled[i];
217 delete fHistw[i];
218 delete fGraphw[i];
219 }
220 delete fHistled;
221 delete fHistw;
222}
223
224void MCaos::ShowHistograms()
225{
226 if (!fHistpr || fHistpr->GetEntries()<1)
227 return;
228
229 TH1 *h;
230
231 TCanvas *c = new TCanvas("cring", "Center of the ring", 800, 800);
232 c->Divide(2,2);
233 c->cd(1);
234 h = (TH1*)fHistprxpry->ProfileX();
235 h->Fit("gaus");
236 h->Draw();
237 h->SetBit(kCanDelete);
238 c->cd(2);
239 h = (TH1*)fHistprxpry->ProfileY();
240 h->Fit("gaus");
241 h->Draw();
242 h->SetBit(kCanDelete);
243 c->cd(3);
244 fHistpr->Fit("gaus");
245 fHistpr->DrawCopy();
246 c->cd(4);
247 fHistprxpry->DrawCopy(/*"surf2"*/);
248 c->Update();
249
250 const Int_t n1 = (Int_t)(sqrt(fPositions.GetEntriesFast())+1.0);
251 const Int_t n2 = (Int_t)(sqrt(fPositions.GetEntriesFast())+0.5);
252
253 TCanvas *c1 = new TCanvas("cpos", "Led Positions", 800, 600);
254 TCanvas *c2 = new TCanvas("cangle", "Angles of the Leds", 800, 600);
255 c1->Divide(n1, n2);
256 c2->Divide(n1, n2);
257 for (int i=0; i<fPositions.GetEntriesFast(); i++)
258 {
259 c1->cd(i+1);
260 fHistled[i]->DrawCopy();
261 c2->cd(i+1);
262 fHistw[i]->DrawCopy();
263 }
264 c1->Update();
265 c2->Update();
266
267 c = new TCanvas("ctime", "timedevelopement of center", 800, 800);
268 c->Divide(1,2);
269 c->cd(1);
270 h = fGraphprx->GetHistogram();
271 h->SetXTitle("time [s]");
272 h->SetYTitle("x [pix]");
273 h->DrawCopy();
274 //fGraphprx->DrawClone("LP*")->SetBit(kCanDelete);
275 c->cd(2);
276 h = fGraphpry->GetHistogram();
277 h->SetXTitle("time [s]");
278 h->SetYTitle("y [pix]");
279 h->DrawCopy();
280 //fGraphpry->DrawClone("LP*")->SetBit(kCanDelete);
281 c->Modified();
282 c->Update();
283
284 /* --------------------------------------------------------
285 * CALCULATE OFFSETS!
286 * --------------------------------------------------------
287 Rings r;
288 r.CalcRings(fPositions);
289
290 const Ring &c = r.GetCenter();
291 */
292}
293
294void MCaos::ResetHistograms()
295{
296 if (!fHistpr)
297 return;
298
299 fHistpr->Reset();
300 fHistprxpry->Reset();
301 fHistallw->Reset();
302 for (int i=0; i<6; i++)
303 {
304 fHistled[i]->Reset();
305 fHistw[i]->Reset();
306 }
307}
308
309void MCaos::Run(byte *img, bool printl, bool printr, const ZdAz &pos, const MTime &t)
310{
311 Leds &leds = *fLeds;
312 leds.Clear();
313
314 // img width height radius sigma
315 FilterLed f(img, 768, 576, 50, 3.0);
316
317 Int_t first=0;
318 for (int i=0; i<fPositions.GetEntriesFast(); i++)
319 {
320 // Try to find Led in this area
321 const Led &l0 = fPositions(i);
322 f.Execute(leds, l0.GetX(), l0.GetY());
323
324 // Loop over newly found Leds
325 for (int j=first; j<leds.GetEntries(); j++)
326 {
327 Led &l1 = leds(j);
328
329 // Add Offset to Led
330 l1.AddOffset(l0.GetDx(), l0.GetDy());
331
332 // Mark Led in image (FIXME: Move to MStarguider)
333 f.MarkPoint(l1.GetX(), l1.GetY(), l1.GetMag());
334
335 // Fill values into Histogram
336 if (!fHistpr)
337 continue;
338
339 fHistled[i]->Fill(l1.GetX(), l1.GetY());
340 fHistw[i]->Fill(l1.GetPhi());
341 }
342 first = leds.GetEntries();
343 }
344
345 Rings rings;
346 rings.CalcRings(leds, 266, 272);
347
348 const Ring &center = rings.GetCenter();
349
350 f.DrawCircle(center, 0x80);
351 f.DrawCircle(center, 5.0, 0x80);
352 f.DrawCircle(center, 115.0, 0x80);
353 f.DrawCircle(center, 190.0, 0x80);
354 // f.MarkPoint(center.GetX(), center.GetY(), 0xa0);
355
356 // FIXME!
357 static const MTime t0(t);
358 fEvtTime = t-t0;
359
360 if (fHistpr)
361 {
362 fHistpr->Fill(center.GetR());
363 fHistprxpry->Fill(center.GetX(), center.GetY());
364 fGraphprx->SetPoint(fGraphprx->GetN(), fEvtTime, center.GetX());
365 fGraphpry->SetPoint(fGraphpry->GetN(), fEvtTime, center.GetY());
366 }
367
368 if (printl)
369 leds.Print();
370 if (printr)
371 rings.Print();
372
373 if (fFile && leds.GetEntries()>0)
374 {
375 fZenithDist = pos.Zd(); //fCosy ? fCosy->GetPointingPos().Zd() : 0
376 fAzimuth = pos.Az(); //fCosy ? fCosy->GetPointingPos().Az() : 0;
377
378 TTree *t = (TTree*)fFile->Get("Data");
379 t->Fill();
380 }
381
382 /*
383 if (fCaosAnalyse->IsEntryEnabled(IDM_kStopAnalyse))
384 {
385 const Ring &center = rings.GetCenter();
386
387 Double_t phi[6] =
388 {
389 -124.727,
390 -61.0495,
391 -16.7907,
392 49.3119,
393 139.086
394 };
395
396 Double_t sum = 0;
397 for (int i=0; i<6 && leds.At(i); i++)
398 {
399 const Double_t w = (leds(i).GetPhi()-phi[i])*60;
400
401 sum += w;
402
403 fHistw[i]->Fill(w);
404 fHistv[i]->Fill(leds(i).GetPhi());
405 fGraphw[i]->SetPoint(fGraphw[i]->GetN(), fEvtTime, w);
406 }
407 fHistallw->Fill(sum/5);
408 }
409 */
410}
411
Note: See TracBrowser for help on using the repository browser.