source: trunk/Cosy/main/MCaos.cc@ 16318

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