Index: trunk/FACT++/drive/MCaos.cc =================================================================== --- trunk/FACT++/drive/MCaos.cc (revision 18621) +++ trunk/FACT++/drive/MCaos.cc (revision 18622) @@ -2,22 +2,9 @@ #include - -#include -#include -#include -#include -#include -#include -#include - -#include "MLog.h" -#include "MLogManip.h" - -#include "MTime.h" -#include "MPointing.h" +#include +#include +#include #include "Led.h" -#include "Ring.h" -#include "Rings.h" #include "FilterLed.h" @@ -29,492 +16,120 @@ if (!fin) { - gLog << err << "ERROR - Cannot open " << name << endl; + cout << "ERROR - Cannot open " << name << endl; return; } - fPositions.Clear(); + fPositions.clear(); - gLog << all << " Reading " << name << ":" << endl; - gLog << inf << "------------------------------" << endl; + cout << " Reading " << name << ":" << endl; + cout << "------------------------------" << endl; while (1) { - Double_t px, py, ox, oy; - fin >> px >> py >> ox >> oy; + double px, py, phi; + fin >> px >> py >> phi; if (!fin) break; - gLog << " Led #" << fPositions.GetEntriesFast() << ": "; - gLog << setw(3) << px << " "; - gLog << setw(3) << py << " ("; - gLog << setw(3) << ox << ", "; - gLog << setw(3) << oy << ")" << endl; - AddPosition(px, py, ox, oy); + cout << " Led #" << fPositions.size() << ": "; + cout << setw(3) << px << " "; + cout << setw(3) << py << " ("; + cout << setw(3) << phi << ")\n"; + AddPosition(px, py, phi); } - gLog << all << "Found " << fPositions.GetEntriesFast() << " leds." << endl; + cout << "Found " << fPositions.size() << " leds." << endl; } -void MCaos::OpenFile() +void MCaos::CalcCenters(const vector &leds, float min, float max) { - int i=0; - char name[100]; - while (1) - { - sprintf(name, "data/data%03d.root", i++); - if (gSystem->AccessPathName(name, kFileExists)) - break; - } + fRings.clear(); - if (fFile) - delete fFile; + const int nPoints = leds.size(); - fFile = new TFile(name, "RECREATE"); + // A minimum of at least 3 points is mandatory! + if (nPointsIsOpen()) - { - delete fFile; - fFile = NULL; +// ofstream fout("rings.txt", ios::app); - gLog << err << "ERROR - Cannot open file '" << name << "'" << endl; - } + for (int i=0; i=0&&ring.GetR()=0&&ring.GetR()>max)) + continue; - tree->Branch("Leds", "TClonesArray", &fLeds); - tree->Branch("ZenithDist.", &fZenithDist, "fZenithDist/D"); - tree->Branch("Azimuth.", &fAzimuth, "fAzimuth/D"); - tree->Branch("EvtTime.", &fEvtTime, "fEvtTime/D"); - - gLog << inf << "Root file '" << name << "' open." << endl; + fRings.push_back(ring); + } } -void MCaos::CloseFile() +int32_t MCaos::CalcRings(std::vector &leds, float min, float max) { - if (!fFile) - return; + CalcCenters(leds, min, max); - const TString name = fFile->GetName(); - const Double_t n = ((TTree*)fFile->Get("Data"))->GetEntries(); + fCenter.InterpolCenters(fRings); - fFile->Write(); - delete fFile; - fFile = NULL; + for (auto it=leds.begin(); it!=leds.end(); it++) + it->CalcPhi(fCenter); - gLog << inf << "Root file closed (n=" << n << ")" << endl; - - if (n<1) - { - gSystem->Unlink(name); - gLog << warn << "Root file deleted - no entries." << endl; - } + return fRings.size(); } -void MCaos::InitHistograms() +Ring MCaos::Run(uint8_t *img) { - if (fHistpr) - return; - - Rings r; - r.SetMinNumberLeds(fMinNumberLeds); - r.CalcRings(fPositions); - - const Ring &c = r.GetCenter(); - - Double_t xmin = c.GetX()-50; - Double_t xmax = c.GetX()+50; - - Double_t ymin = c.GetY()-50; - Double_t ymax = c.GetY()+50; - - Double_t rmin = c.GetR()-50; - Double_t rmax = c.GetR()+50; - - Int_t xbin = 1001; - Int_t ybin = 1001; - Int_t rbin = 1001; - - const Int_t sz = 50; - fHistled = new TH2F*[fPositions.GetEntriesFast()]; - fHistw = new TH1F*[fPositions.GetEntriesFast()]; - for (int i=0; iSetXTitle("x [pix]"); - fHistled[i]->SetYTitle("counts"); - - name = "Angle"; - title = "Angle of the Led #"; - - name += i; - title += i; - - fHistw[i] = new TH1F; - fHistw[i]->SetNameTitle(name, title); - fHistw[i]->SetBins(101, -50.5, 50.5); - fHistw[i]->SetXTitle("\\Phi [arcmin]"); - fHistw[i]->SetYTitle("counts"); - } - - fHistallw = new TH1F; - fHistallw->SetNameTitle("allw","Rotation angel"); - fHistallw->SetBins(26, -25, 25); - fHistallw->SetXTitle("\\phi [arcmin]"); - fHistallw->SetYTitle("counts"); - - fHistprxpry = new TH2F; - fHistprxpry->SetNameTitle("prx und pry","x- and y-coordinate of the ring-center"); - fHistprxpry->SetBins(xbin, xmin, xmax, ybin, ymin, ymax); - fHistprxpry->SetXTitle("x [pix]"); - fHistprxpry->SetYTitle("y [pix]"); - fHistprxpry->SetZTitle("counts"); - - fGraphprx = new TGraph; - fGraphprx->SetTitle("time-developement of the x-coordinate of the ring-center"); - - fGraphpry = new TGraph; - fGraphpry->SetTitle("time-developement of the y-coordinate of the ring-center"); - - fGraphw = new TGraph; - fGraphw->SetTitle("Time-developement of rotation angle"); - - fHistpr = new TH1F("pr","Radius of the ring", rbin, rmin, rmax); - fHistpr->SetXTitle("r [pix]"); - fHistpr->SetYTitle("counts"); -} - -void MCaos::DeleteHistograms() -{ - TH1F *dummy = fHistpr; - fHistpr=NULL; - - if (!dummy) - return; - - delete dummy; - delete fHistprxpry; - delete fHistallw; - delete fGraphprx; - delete fGraphpry; - delete fGraphr; - - for (int i=0; i<6; i++) - { - delete fHistled[i]; - delete fHistw[i]; - delete fGraphw; - } - delete fHistled; - delete fHistw; -} - -void MCaos::ShowHistograms() -{ - if (!fHistpr || fHistpr->GetEntries()<1) - return; - - TH1 *h; - - TCanvas *c = new TCanvas("cring", "Center of the ring", 800, 800); - c->Divide(2,2); - c->cd(1); - h = (TH1*)fHistprxpry->ProfileX(); - h->Fit("gaus"); - h->Draw(); - h->SetBit(kCanDelete); - c->cd(2); - h = (TH1*)fHistprxpry->ProfileY(); - h->Fit("gaus"); - h->Draw(); - h->SetBit(kCanDelete); - c->cd(3); - fHistpr->Fit("gaus"); - fHistpr->DrawCopy(); - c->cd(4); - fHistprxpry->DrawCopy(/*"surf2"*/); - c->Update(); - - const Int_t n1 = (Int_t)(TMath::Sqrt(fPositions.GetEntriesFast())+1.0); - const Int_t n2 = (Int_t)(TMath::Sqrt(fPositions.GetEntriesFast())+0.5); - - TCanvas *c1 = new TCanvas("cpos", "Led Positions", 800, 600); - TCanvas *c2 = new TCanvas("cangle", "Angles of the Leds", 800, 600); - c1->Divide(n1, n2); - c2->Divide(n1, n2); - for (int i=0; icd(i+1); - fHistled[i]->DrawCopy(); - c2->cd(i+1); - fHistw[i]->DrawCopy(); - } - c1->Update(); - c2->Update(); - - /* - c = new TCanvas("ctime", "Timedevelopement of Center", 800, 800); - c->Divide(1,3); - c->cd(1); - h = fGraphprx->GetHistogram(); - h->SetXTitle("time [s]"); - h->SetYTitle("x [pix]"); - h->DrawCopy(); - c->SetSelectedPad(NULL); - fGraphprx->DrawClone("ALP*")->SetBit(kCanDelete); - gPad->Modified(); - gPad->Update(); - c->cd(2); - h = fGraphpry->GetHistogram(); - h->SetXTitle("time [s]"); - h->SetYTitle("y [pix]"); - h->DrawCopy(); - //((TPad*)gPad)->SetSelectedPad(NULL); - //fGraphpry->DrawClone("ALP*")->SetBit(kCanDelete); - c->cd(3); - h = fGraphr->GetHistogram(); - h->SetXTitle("time [s]"); - h->SetYTitle("r [pix]"); - h->DrawCopy(); - //((TPad*)gPad)->SetSelectedPad(NULL); - //fGraphr->DrawClone("ALP*")->SetBit(kCanDelete); - c->Modified(); - c->Update(); - */ - - c = new TCanvas("crot", "rotation angle", 800, 600); - c->Divide(2,1); - c->cd(1); - fHistallw->SetXTitle("\\phi [arcmin]"); - fHistallw->SetYTitle("counts"); - fHistallw->DrawCopy(); - /* - c->cd(2); - h = fGraphw->GetHistogram(); - h->SetXTitle("time [s]"); - h->SetYTitle("\\phi [arcmin]"); - h->DrawCopy(); - ((TPad*)gPad)->SetSelected(NULL); - fGraphw->DrawClone("ALP*")->SetBit(kCanDelete); - */ - - /* -------------------------------------------------------- - * CALCULATE OFFSETS! - * -------------------------------------------------------- - Rings r; - r.CalcRings(fPositions); - - const Ring &c = r.GetCenter(); - */ -} - -void MCaos::ResetHistograms() -{ - if (!fHistpr) - return; - - fHistpr->Reset(); - fHistprxpry->Reset(); - fHistallw->Reset(); - for (int i=0; i<6; i++) - { - fHistled[i]->Reset(); - fHistw[i]->Reset(); - } -} - -Ring MCaos::Run(byte *img, bool printl, bool printr, const ZdAz &pos, const MTime &t) -{ - Leds &leds = *fLeds; - leds.Clear(); - - fNumDetectedLEDs = 0; - fNumDetectedRings = 0; - - /* - //the following lines are just to test the new setup - static int i=0; - i++; - i%=2; - if (i==0) - ReadResources("leds0.txt"); - else - ReadResources("leds1.txt"); - cout << "LEDs " << i << " " << flush; - */ + fLeds.clear(); // img width height radius sigma FilterLed f(img, 768, 576, fSizeBox, fSizeBox, fCut); - Int_t first=0; - for (int i=0; i arr; + // Try to find Led in this area - const Led &l0 = fPositions(i); - f.Execute(leds, TMath::FloorNint(l0.GetX()), TMath::FloorNint(l0.GetY())); - - fNumDetectedLEDs = leds.GetEntries(); + f.Execute(arr, floor(it->GetX()), floor(it->GetY())); // Loop over newly found Leds - for (int j=first; jAddOffset(it->GetDx(), it->GetDy()); - // Add Offset to Led - l1.AddOffset(l0.GetDx(), l0.GetDy()); + // Remember the expected phi for each detected led + jt->SetPhi(it->GetPhi()); // Mark Led in image (FIXME: Move to MStarguider) - f.MarkPoint(l1.GetX(), l1.GetY(), l1.GetMag()); + f.MarkPoint(jt->GetX(), jt->GetY(), jt->GetMag()); + } - //old - /* - // Fill values into Histogram - if (!fHistpr) - continue; - - fHistled[i]->Fill(l1.GetX(), l1.GetY()); - fHistw[i]->Fill(l1.GetPhi()); - */ - } - first = leds.GetEntries(); + fLeds.insert(fLeds.end(), arr.begin(), arr.end()); } - Rings rings; - rings.SetMinNumberLeds(fMinNumberLeds); -// rings.CalcRings(leds, 265, 268); -// rwagner -// rings.CalcRings(leds, 158, 164); - fNumDetectedRings = rings.CalcRings(leds, fMinRadius, fMaxRadius); + fNumDetectedRings = CalcRings(fLeds, fMinRadius, fMaxRadius); - const Ring ¢er = rings.GetCenter(); + double sumphi = 0; + for (auto it=fLeds.begin(); it!=fLeds.end(); it++) + { + //cout << it->CalcPhi(fCenter) << "|" << it->GetPhi() << " "; + double dphi = it->CalcPhi(fCenter) - it->GetPhi(); + if (dphi>M_PI) + dphi -= 2*M_PI; + if (dphi<-M_PI) + dphi += 2*M_PI; -//uncommented for testing -// center.Print(); + sumphi += dphi; + } + //cout << endl; - // FIXME! - static const MTime t0(t); - fEvtTime = t-t0; + fCenter.SetPhi(sumphi/fLeds.size()); - if (fHistpr) - { - fHistpr->Fill(center.GetR()); - fHistprxpry->Fill(center.GetX(), center.GetY()); - fGraphprx->SetPoint(fGraphprx->GetN(), fEvtTime, center.GetX()); - fGraphpry->SetPoint(fGraphpry->GetN(), fEvtTime, center.GetY()); - - //new - //----- - Double_t sum = 0; - for (int j=0; jFill(l1.GetX(), l1.GetY()); - //fHistw[j]->Fill(l1.GetPhi()); - - Double_t phi[6] = - { - 0, - 0, - 0, - 0, - 0, - 0 - }; - - const Double_t w = (leds(j).GetPhi()-phi[j])*60; - sum += w; - - fHistw[j]->Fill(w); - sum /= leds.GetEntries(); - } - fGraphw->SetPoint(fGraphw->GetN(), fEvtTime, sum); - fHistallw->Fill(sum/leds.GetEntries()); - //----- - } - - /* - //test - give number of rings - cout << rings.GetEntries() << " " << flush; - */ - - if (printl) - leds.Print(); - if (printr) - rings.Print(); - - if (fFile && leds.GetEntries()>0) - { - fZenithDist = pos.Zd(); //fCosy ? fCosy->GetPointingPos().Zd() : 0 - fAzimuth = pos.Az(); //fCosy ? fCosy->GetPointingPos().Az() : 0; - - TTree *t = (TTree*)fFile->Get("Data"); - t->Fill(); - } - - return center; - /* - if (fCaosAnalyse->IsEntryEnabled(IDM_kStopAnalyse)) - { - const Ring ¢er = rings.GetCenter(); - - Double_t phi[6] = - { - -124.727, - -61.0495, - -16.7907, - 49.3119, - 139.086 - }; - - Double_t sum = 0; - for (int i=0; i<6 && leds.At(i); i++) - { - const Double_t w = (leds(i).GetPhi()-phi[i])*60; - - sum += w; - - fHistw[i]->Fill(w); - fHistv[i]->Fill(leds(i).GetPhi()); - fGraphw[i]->SetPoint(fGraphw[i]->GetN(), fEvtTime, w); - } - fHistallw->Fill(sum/5); - } - */ + return fCenter; } - - -Int_t MCaos::ReadEnv(const TEnv &env, TString prefix, Bool_t print) -{ - if (IsEnvDefined(env, prefix, "File", print)) - ReadResources(GetEnvValue(env, prefix, "File", "ledsxxx.txt")); - - if (IsEnvDefined(env, prefix, "RadiusMin", print)) - fMinRadius = GetEnvValue(env, prefix, "RadiusMin", fMinRadius); - - if (IsEnvDefined(env, prefix, "RadiusMax", print)) - fMaxRadius = GetEnvValue(env, prefix, "RadiusMax", fMaxRadius); - - if (IsEnvDefined(env, prefix, "MinNumberLeds", print)) - fMinNumberLeds = GetEnvValue(env, prefix, "MinNumberLeds", fMinNumberLeds); - - if (IsEnvDefined(env, prefix, "SizeBox", print)) - fSizeBox = GetEnvValue(env, prefix, "SizeBox", fSizeBox); - - if (IsEnvDefined(env, prefix, "CleaningLevel", print)) - fCut = GetEnvValue(env, prefix, "CleaningLevel", fCut); - - if (IsEnvDefined(env, prefix, "ArcsecPerPixel", print)) - fArcsecPerPixel = GetEnvValue(env, prefix, "ArcsecPerPixel", fArcsecPerPixel); - - return kTRUE; -} Index: trunk/FACT++/drive/MCaos.h =================================================================== --- trunk/FACT++/drive/MCaos.h (revision 18621) +++ trunk/FACT++/drive/MCaos.h (revision 18622) @@ -2,107 +2,57 @@ #define CAOS_MCaos -#ifndef MARS_MParContainer -#include "MParContainer.h" -#endif -#ifndef CAOS_Leds -#include "Leds.h" -#endif #ifndef CAOS_Ring #include "Ring.h" #endif -#ifndef COSY_PixClient -#include "PixClient.h" // byte -#endif - -class TFile; -class TH1F; -class TH2F; -class TGraph; - -class MTime; -class ZdAz; - -class MCaos : public MParContainer +class MCaos { private: - TFile *fFile; // File we may write data to + std::vector fPositions; + std::vector fLeds; - Leds fPositions; + int16_t fMinNumberLeds; // minimum number of detected leds required + double fMinRadius; // minimum radius for cut in ring radius + double fMaxRadius; // maximum radius for cut in ring radius + uint16_t fSizeBox; // Size of the search box (side length in units of pixels) + double fCut; // Cleaning level (sigma above noise) - Leds *fLeds; - Double_t fEvtTime; - Double_t fZenithDist; - Double_t fAzimuth; + int32_t fNumDetectedRings; - TH1F *fHistpr; - TH2F **fHistled; - TH1F *fHistallw; - TH1F **fHistw; + Ring fCenter; + std::vector fRings; - TH2F *fHistprxpry; - - TGraph *fGraphprx; - TGraph *fGraphpry; - TGraph *fGraphw; - TGraph *fGraphr; - - Short_t fMinNumberLeds; // minimum number of detected leds required - Double_t fMinRadius; // minimum radius for cut in ring radius - Double_t fMaxRadius; // maximum radius for cut in ring radius - UShort_t fSizeBox; // Size of the search box (side length in units of pixels) - Double_t fCut; // Cleaning level (sigma above noise) - Double_t fArcsecPerPixel; // Conversion from arcseconds to pixel - - Int_t fNumDetectedLEDs; - Int_t fNumDetectedRings; + void CalcCenters(const std::vector &leds, float min, float max); + int32_t CalcRings(std::vector &leds, float min=-1, float max=-1); + const Ring &GetCenter() const { return fCenter; } public: - MCaos(const char *name=0, const char *title=0) - : fFile(NULL), fHistpr(NULL), fMinNumberLeds(5), - fMinRadius(265), fMaxRadius(268), fSizeBox(19), fCut(3.0) + MCaos() : fMinRadius(236.7), fMaxRadius(238.6), fSizeBox(19), fCut(3.5) { - fLeds = new Leds; } ~MCaos() { - CloseFile(); - DeleteHistograms(); - delete fLeds; } - void AddPosition(Float_t x, Float_t y, Float_t dx, Float_t dy) + void AddPosition(float x, float y, float phi) { - fPositions.Add(x, y, dx, dy); + fPositions.push_back(Led(x, y, phi)); } void ReadResources(const char *name="leds.txt"); - void OpenFile(); - void CloseFile(); - - void SetMinNumberLeds(Short_t n) + void SetMinNumberLeds(int16_t n) { fMinNumberLeds = n; } - void SetMinRadius(Double_t min) { fMinRadius=min; } - void SetMaxRadius(Double_t max) { fMaxRadius=max; } + void SetMinRadius(double min) { fMinRadius=min; } + void SetMaxRadius(double max) { fMaxRadius=max; } - void InitHistograms(); - void DeleteHistograms(); - void ShowHistograms(); - void ResetHistograms(); + int32_t GetNumDetectedLEDs() const { return fLeds.size(); } + int32_t GetNumDetectedRings() const { return fNumDetectedRings; } - Int_t GetNumDetectedLEDs() const { return fNumDetectedLEDs; } - Int_t GetNumDetectedRings() const { return fNumDetectedRings; } - - Double_t GetArcsecPerPixel() const { return fArcsecPerPixel; } - - Ring Run(byte *img, bool printl, bool printr, const ZdAz &pos, - const MTime &t); - - Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print=kFALSE); + Ring Run(uint8_t *img); };