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

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