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

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