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

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