#include "MCaos.h" #include #include #include #include #include #include #include #include #include #include "MTime.h" #include "Led.h" #include "Ring.h" #include "Rings.h" #include "FilterLed.h" #include "coord.h" using namespace std; void MCaos::ReadResources(const char *name) { ifstream fin(name); if (!fin) { cout << "ERROR - Cannot open " << name << endl; return; } fPositions.Clear(); cout << " Reading " << name << ":" << endl; cout << "------------------------------" << endl; while (1) { Double_t px, py, ox, oy; fin >> px >> py >> ox >> oy; if (!fin) break; cout << " Led #" << fPositions.GetEntriesFast() << ": "; cout << setw(3) << px << " "; cout << setw(3) << py << " ("; cout << setw(3) << ox << ", "; cout << setw(3) << oy << ")" << endl; AddPosition(px, py, ox, oy); } cout << "Found " << fPositions.GetEntriesFast() << " leds." << endl; } void MCaos::OpenFile() { int i=0; char name[100]; while (1) { sprintf(name, "data/data%03d.root", i++); if (gSystem->AccessPathName(name, kFileExists)) break; } if (fFile) delete fFile; fFile = new TFile(name, "RECREATE"); if (!fFile->IsOpen()) { delete fFile; fFile = NULL; cout << "Error: Cannot open file '" << name << "'" << endl; } TTree *tree = new TTree("Data", "Real CaOs Data"); fLeds = new Leds; fEvtTime = 0; tree->Branch("Leds", "TClonesArray", &fLeds); tree->Branch("ZenithDist.", &fZenithDist, "fZenithDist/D"); tree->Branch("Azimuth.", &fAzimuth, "fAzimuth/D"); tree->Branch("EvtTime.", &fEvtTime, "fEvtTime/D"); cout << "Root file '" << name << "' open." << endl; } void MCaos::CloseFile() { if (!fFile) return; const TString name = fFile->GetName(); const Double_t n = ((TTree*)fFile->Get("Data"))->GetEntries(); fFile->Write(); delete fFile; fFile = NULL; cout << "Root file closed (n=" << n << ")" << endl; if (n<1) { gSystem->Unlink(name); cout << "Root file deleted - no entries." << endl; } } void MCaos::InitHistograms() { 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, Int_t box, Double_t cut) { Leds &leds = *fLeds; leds.Clear(); /* //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; */ // img width height radius sigma FilterLed f(img, 768, 576, box, cut); Int_t first=0; for (int i=0; iFill(l1.GetX(), l1.GetY()); fHistw[i]->Fill(l1.GetPhi()); */ } first = leds.GetEntries(); } Rings rings; rings.SetMinNumberLeds(fMinNumberLeds); // rings.CalcRings(leds, 265, 268); // rwagner // rings.CalcRings(leds, 158, 164); rings.CalcRings(leds, fMinRadius, fMaxRadius); const Ring ¢er = rings.GetCenter(); //uncommented for testing // center.Print(); // FIXME! static const MTime t0(t); fEvtTime = t-t0; 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); } */ }