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

Last change on this file since 8015 was 7790, checked in by tbretz, 19 years ago
*** empty log message ***
File size: 11.5 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 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.SetMinNumberLeds(fMinNumberLeds);
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(101, -50.5, 50.5);
162 fHistw[i]->SetXTitle("\\Phi [arcmin]");
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-coordinate of the ring-center");
174 fHistprxpry->SetBins(xbin, xmin, xmax, ybin, ymin, ymax);
175 fHistprxpry->SetXTitle("x [pix]");
176 fHistprxpry->SetYTitle("y [pix]");
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;
186 fGraphw->SetTitle("Time-developement of rotation angle");
187
188 fHistpr = new TH1F("pr","Radius of the ring", rbin, rmin, rmax);
189 fHistpr->SetXTitle("r [pix]");
190 fHistpr->SetYTitle("counts");
191}
192
193void MCaos::DeleteHistograms()
194{
195 TH1F *dummy = fHistpr;
196 fHistpr=NULL;
197
198 if (!dummy)
199 return;
200
201 delete dummy;
202 delete fHistprxpry;
203 delete fHistallw;
204 delete fGraphprx;
205 delete fGraphpry;
206 delete fGraphr;
207
208 for (int i=0; i<6; i++)
209 {
210 delete fHistled[i];
211 delete fHistw[i];
212 delete fGraphw;
213 }
214 delete fHistled;
215 delete fHistw;
216}
217
218void MCaos::ShowHistograms()
219{
220 if (!fHistpr || fHistpr->GetEntries()<1)
221 return;
222
223 TH1 *h;
224
225 TCanvas *c = new TCanvas("cring", "Center of the ring", 800, 800);
226 c->Divide(2,2);
227 c->cd(1);
228 h = (TH1*)fHistprxpry->ProfileX();
229 h->Fit("gaus");
230 h->Draw();
231 h->SetBit(kCanDelete);
232 c->cd(2);
233 h = (TH1*)fHistprxpry->ProfileY();
234 h->Fit("gaus");
235 h->Draw();
236 h->SetBit(kCanDelete);
237 c->cd(3);
238 fHistpr->Fit("gaus");
239 fHistpr->DrawCopy();
240 c->cd(4);
241 fHistprxpry->DrawCopy(/*"surf2"*/);
242 c->Update();
243
244 const Int_t n1 = (Int_t)(sqrt(fPositions.GetEntriesFast())+1.0);
245 const Int_t n2 = (Int_t)(sqrt(fPositions.GetEntriesFast())+0.5);
246
247 TCanvas *c1 = new TCanvas("cpos", "Led Positions", 800, 600);
248 TCanvas *c2 = new TCanvas("cangle", "Angles of the Leds", 800, 600);
249 c1->Divide(n1, n2);
250 c2->Divide(n1, n2);
251 for (int i=0; i<fPositions.GetEntriesFast(); i++)
252 {
253 c1->cd(i+1);
254 fHistled[i]->DrawCopy();
255 c2->cd(i+1);
256 fHistw[i]->DrawCopy();
257 }
258 c1->Update();
259 c2->Update();
260
261 /*
262 c = new TCanvas("ctime", "Timedevelopement of Center", 800, 800);
263 c->Divide(1,3);
264 c->cd(1);
265 h = fGraphprx->GetHistogram();
266 h->SetXTitle("time [s]");
267 h->SetYTitle("x [pix]");
268 h->DrawCopy();
269 c->SetSelectedPad(NULL);
270 fGraphprx->DrawClone("ALP*")->SetBit(kCanDelete);
271 gPad->Modified();
272 gPad->Update();
273 c->cd(2);
274 h = fGraphpry->GetHistogram();
275 h->SetXTitle("time [s]");
276 h->SetYTitle("y [pix]");
277 h->DrawCopy();
278 //((TPad*)gPad)->SetSelectedPad(NULL);
279 //fGraphpry->DrawClone("ALP*")->SetBit(kCanDelete);
280 c->cd(3);
281 h = fGraphr->GetHistogram();
282 h->SetXTitle("time [s]");
283 h->SetYTitle("r [pix]");
284 h->DrawCopy();
285 //((TPad*)gPad)->SetSelectedPad(NULL);
286 //fGraphr->DrawClone("ALP*")->SetBit(kCanDelete);
287 c->Modified();
288 c->Update();
289 */
290
291 c = new TCanvas("crot", "rotation angle", 800, 600);
292 c->Divide(2,1);
293 c->cd(1);
294 fHistallw->SetXTitle("\\phi [arcmin]");
295 fHistallw->SetYTitle("counts");
296 fHistallw->DrawCopy();
297 /*
298 c->cd(2);
299 h = fGraphw->GetHistogram();
300 h->SetXTitle("time [s]");
301 h->SetYTitle("\\phi [arcmin]");
302 h->DrawCopy();
303 ((TPad*)gPad)->SetSelected(NULL);
304 fGraphw->DrawClone("ALP*")->SetBit(kCanDelete);
305 */
306
307 /* --------------------------------------------------------
308 * CALCULATE OFFSETS!
309 * --------------------------------------------------------
310 Rings r;
311 r.CalcRings(fPositions);
312
313 const Ring &c = r.GetCenter();
314 */
315}
316
317void MCaos::ResetHistograms()
318{
319 if (!fHistpr)
320 return;
321
322 fHistpr->Reset();
323 fHistprxpry->Reset();
324 fHistallw->Reset();
325 for (int i=0; i<6; i++)
326 {
327 fHistled[i]->Reset();
328 fHistw[i]->Reset();
329 }
330}
331
332Ring MCaos::Run(byte *img, bool printl, bool printr, const ZdAz &pos, const MTime &t, Double_t box, Double_t cut)
333{
334 Leds &leds = *fLeds;
335 leds.Clear();
336
337 /*
338 //the following lines are just to test the new setup
339 static int i=0;
340 i++;
341 i%=2;
342 if (i==0)
343 ReadResources("leds0.txt");
344 else
345 ReadResources("leds1.txt");
346 cout << "LEDs " << i << " " << flush;
347 */
348
349 // img width height radius sigma
350 FilterLed f(img, 768, 576, box, cut);
351
352 Int_t first=0;
353 for (int i=0; i<fPositions.GetEntriesFast(); i++)
354 {
355 // Try to find Led in this area
356 const Led &l0 = fPositions(i);
357 f.Execute(leds, l0.GetX(), l0.GetY());
358
359 // Loop over newly found Leds
360 for (int j=first; j<leds.GetEntries(); j++)
361 {
362 Led &l1 = leds(j);
363
364 // Add Offset to Led
365 l1.AddOffset(l0.GetDx(), l0.GetDy());
366
367 // Mark Led in image (FIXME: Move to MStarguider)
368 f.MarkPoint(l1.GetX(), l1.GetY(), l1.GetMag());
369
370 //old
371 /*
372 // Fill values into Histogram
373 if (!fHistpr)
374 continue;
375
376 fHistled[i]->Fill(l1.GetX(), l1.GetY());
377 fHistw[i]->Fill(l1.GetPhi());
378 */
379 }
380 first = leds.GetEntries();
381 }
382
383 Rings rings;
384 rings.SetMinNumberLeds(fMinNumberLeds);
385// rings.CalcRings(leds, 265, 268);
386// rwagner
387// rings.CalcRings(leds, 158, 164);
388 rings.CalcRings(leds, fMinRadius, fMaxRadius);
389
390 const Ring &center = rings.GetCenter();
391
392//uncommented for testing
393// center.Print();
394
395 // FIXME!
396 static const MTime t0(t);
397 fEvtTime = t-t0;
398
399 if (fHistpr)
400 {
401 fHistpr->Fill(center.GetR());
402 fHistprxpry->Fill(center.GetX(), center.GetY());
403 fGraphprx->SetPoint(fGraphprx->GetN(), fEvtTime, center.GetX());
404 fGraphpry->SetPoint(fGraphpry->GetN(), fEvtTime, center.GetY());
405
406 //new
407 //-----
408 Double_t sum = 0;
409 for (int j=0; j<leds.GetEntries(); j++)
410 {
411 Led &l1 = leds(j);
412
413 fHistled[j]->Fill(l1.GetX(), l1.GetY());
414 //fHistw[j]->Fill(l1.GetPhi());
415
416 Double_t phi[6] =
417 {
418 0,
419 0,
420 0,
421 0,
422 0,
423 0
424 };
425
426 const Double_t w = (leds(j).GetPhi()-phi[j])*60;
427 sum += w;
428
429 fHistw[j]->Fill(w);
430 sum /= leds.GetEntries();
431 }
432 fGraphw->SetPoint(fGraphw->GetN(), fEvtTime, sum);
433 fHistallw->Fill(sum/leds.GetEntries());
434 //-----
435 }
436
437 /*
438 //test - give number of rings
439 cout << rings.GetEntries() << " " << flush;
440 */
441
442 if (printl)
443 leds.Print();
444 if (printr)
445 rings.Print();
446
447 if (fFile && leds.GetEntries()>0)
448 {
449 fZenithDist = pos.Zd(); //fCosy ? fCosy->GetPointingPos().Zd() : 0
450 fAzimuth = pos.Az(); //fCosy ? fCosy->GetPointingPos().Az() : 0;
451
452 TTree *t = (TTree*)fFile->Get("Data");
453 t->Fill();
454 }
455
456 return center;
457 /*
458 if (fCaosAnalyse->IsEntryEnabled(IDM_kStopAnalyse))
459 {
460 const Ring &center = rings.GetCenter();
461
462 Double_t phi[6] =
463 {
464 -124.727,
465 -61.0495,
466 -16.7907,
467 49.3119,
468 139.086
469 };
470
471 Double_t sum = 0;
472 for (int i=0; i<6 && leds.At(i); i++)
473 {
474 const Double_t w = (leds(i).GetPhi()-phi[i])*60;
475
476 sum += w;
477
478 fHistw[i]->Fill(w);
479 fHistv[i]->Fill(leds(i).GetPhi());
480 fGraphw[i]->SetPoint(fGraphw[i]->GetN(), fEvtTime, w);
481 }
482 fHistallw->Fill(sum/5);
483 }
484 */
485}
486
Note: See TracBrowser for help on using the repository browser.