#include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "coord.h" #include "MGList.h" #include "MBending.h" class Set : public TObject { friend ifstream &operator>>(ifstream &fin, Set &set); private: Double_t fStarAz; Double_t fStarEl; Double_t fRawAz; Double_t fRawEl; public: Set(Double_t sel=0, Double_t saz=0, Double_t rel=0, Double_t raz=0) : fStarAz(saz*TMath::Pi()/180), fStarEl(sel*TMath::Pi()/180), fRawAz(raz*TMath::Pi()/180), fRawEl(rel*TMath::Pi()/180) { } Double_t GetResidual() const { Double_t del = fRawEl-fStarEl; Double_t daz = fRawAz-fStarAz; Double_t dphi2 = daz/2.; Double_t cos2 = cos(dphi2)*cos(dphi2); Double_t sin2 = sin(dphi2)*sin(dphi2); Double_t d = cos(del)*cos2 - cos(fRawEl+fStarEl)*sin2; Double_t dist = acos(d); return dist * 180 / TMath::Pi(); } void operator=(Set &set) { fStarAz = set.fStarAz; fStarEl = set.fStarEl; fRawAz = set.fRawAz; fRawEl = set.fRawEl; } Double_t GetDEl() const { return (fRawEl-fStarEl)*180/TMath::Pi(); } Double_t GetDZd() const { return -GetDEl(); } Double_t GetDAz() const { return (fRawAz-fStarAz)*180/TMath::Pi(); } Double_t GetStarEl() const { return fStarEl*180/TMath::Pi(); } Double_t GetStarZd() const { return 90.-fStarEl*180/TMath::Pi(); } Double_t GetStarAz() const { return fStarAz*180/TMath::Pi(); } Double_t GetRawEl() const { return fRawEl*180/TMath::Pi(); } Double_t GetRawAz() const { return fRawAz*180/TMath::Pi(); } Double_t GetRawZd() const { return 90.-fRawEl*180/TMath::Pi(); } ZdAz GetStarZdAz() const { return ZdAz(TMath::Pi()/2-fStarEl, fStarAz); } AltAz GetStarAltAz() const { return AltAz(fStarEl, fStarAz); } ZdAz GetRawZdAz() const { return ZdAz(TMath::Pi()/2-fRawEl, fRawAz); } AltAz GetRawAltAz() const { return AltAz(fRawEl, fRawAz); } void AdjustEl(Double_t del) { fStarEl += del*TMath::Pi()/180; } void AdjustAz(Double_t daz) { fStarAz += daz*TMath::Pi()/180; } void Adjust(const MBending &bend) { AltAz p = bend(GetStarAltAz()); fStarEl = p.Alt(); fStarAz = p.Az(); } void AdjustBack(const MBending &bend) { AltAz p = bend.CorrectBack(GetRawAltAz()); fRawEl = p.Alt(); fRawAz = p.Az(); } ClassDef(Set, 0) }; ClassImp(Set); ifstream &operator>>(ifstream &fin, Set &set) { Double_t v[4]; fin >> v[0]; fin >> v[1]; fin >> v[2]; fin >> v[3]; Double_t dummy; fin>>dummy; fin>>dummy; fin>>dummy; set.fStarAz = v[0]*TMath::Pi()/180; set.fStarEl = v[1]*TMath::Pi()/180; set.fRawAz = v[2]*TMath::Pi()/180; set.fRawEl = v[3]*TMath::Pi()/180; return fin; } class MFit : public TGMainFrame { private: enum { kTbFit = 19, //MBending::GetNumPar(), // FIXME!!! kTbLoad, kTbSave, kTbLoadStars, kTbReset, kTbResetStars }; MGList *fList; TList fOriginal; TList fCoordinates; TList fLabel; MBending fBending; FontStruct_t fFont; void Fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) { f = 0; MBending bend; bend.SetParameters(par); // Set Parameters [deg] to MBending for (int i=0; iGetObjectFit())->Fcn(npar, gin, f, par, iflag); } void DrawMarker(TVirtualPad *pad, Double_t r0, Double_t phi0, Double_t r1, Double_t phi1) { TView *view = pad->GetView(); if (!view) { cout << "No View!" << endl; return; } TMarker mark0; TMarker mark1; mark0.SetMarkerStyle(kStar); mark1.SetMarkerStyle(kStar); mark1.SetMarkerColor(kRed); r0 /= 90; r1 /= 90; phi0 *= TMath::Pi()/180; phi1 *= TMath::Pi()/180; Double_t x0[3] = { r0*cos(phi0), r0*sin(phi0), 0}; Double_t x1[3] = { r1*cos(phi1), r1*sin(phi1), 0}; Double_t y0[3], y1[3]; view->WCtoNDC(x0, y0); view->WCtoNDC(x1, y1); mark0.DrawMarker(y0[0], y0[1]); mark1.DrawMarker(y1[0], y1[1]); } void DrawPolLine(TVirtualPad *pad, Double_t r0, Double_t phi0, Double_t r1, Double_t phi1) { TView *view = pad->GetView(); if (!view) { cout << "No View!" << endl; return; } /* if (r0<0) { r0 = -r0; phi0 += 180; } if (r1<0) { r1 = -r1; phi1 += 180; } phi0 = fmod(phi0+360, 360); phi1 = fmod(phi1+360, 360); if (phi1-phi0<-180) phi1+=360; */ TLine line; line.SetLineWidth(2); line.SetLineColor(kBlue); Double_t p0 = phi0phi1) { Double_t d = r1; r1 = r0; r0 = d; } r0 /= 90; r1 /= 90; Double_t dr = r1-r0; Double_t dp = p1-p0; Double_t x0[3] = { r0*cos(p0*TMath::Pi()/180), r0*sin(p0*TMath::Pi()/180), 0}; for (double i=p0+10; ip1) i=p1; Double_t r = dr/dp*(i-p0)+r0; Double_t p = TMath::Pi()/180*i; Double_t x1[3] = { r*cos(p), r*sin(p), 0}; Double_t y0[3], y1[3]; view->WCtoNDC(x0, y0); view->WCtoNDC(x1, y1); line.DrawLine(y0[0], y0[1], y1[0], y1[1]); x0[0] = x1[0]; x0[1] = x1[1]; } } void DrawSet(TVirtualPad *pad, Set &set, Float_t scale=-1) { Double_t r0 = set.GetRawZd(); Double_t phi0 = set.GetRawAz(); Double_t r1 = set.GetStarZd(); Double_t phi1 = set.GetStarAz(); if (r0<0) { r0 = -r0; phi0 += 180; } if (r1<0) { r1 = -r1; phi1 += 180; } phi0 = fmod(phi0+360, 360); phi1 = fmod(phi1+360, 360); if (phi1-phi0<-180) phi1+=360; if (scale<0 || scale>1000) scale = -1; if (scale>0) { Double_t d = r1-r0; r0 += scale*d; r1 -= scale*d; d = phi1-phi0; phi0 += scale*d; phi1 -= scale*d; DrawPolLine(pad, r0, phi0, r1, phi1); DrawMarker(pad, r0, phi0, r1, phi1); } else DrawMarker(pad, r1, phi1, 0 ,0); } void Fit() { if (fOriginal.GetSize()==0) { cout << "Sorry, no input data loaded..." << endl; return; } fCoordinates.Delete(); for (int i=0; iFindWidget(i); minuit.FixParameter(i); if (l->GetState()==kButtonDown) minuit.Release(i); } //minuit.Command("SHOW PARAMETERS"); //minuit.Command("SHOW LIMITS"); cout << endl; cout << "Starting fit..." << endl; cout << "For the fit an measurement error in the residual of "; cout << "0.02deg (=1SE) is assumed." << endl; cout << endl; Int_t ierflg = 0; ierflg = minuit.Migrad(); cout << "Migrad returns " << ierflg << endl; // minuit.Release(2); ierflg = minuit.Migrad(); cout << "Migrad returns " << ierflg << endl << endl; // // Get Fit Results // fBending.GetMinuitParameters(minuit); fBending.PrintMinuitParameters(minuit); //fBending.Save("bending_magic.txt"); // // Make a copy of all list entries // TList list; list.SetOwner(); for (int i=0; i180) dz -= 360; gdzd.SetPoint(i, za.Az(), set0.GetDZd()); gdaz.SetPoint(i, za.Zd(), dz); graz.SetPoint(i, za.Az(), set0.GetResidual()); grzd.SetPoint(i, za.Zd(), set0.GetResidual()); gaz.SetPoint( i, za.Az(), dz); gzd.SetPoint( i, za.Zd(), set0.GetDZd()); } // // Check for overflows // const Stat_t ov = hres2.GetBinContent(hres2.GetNbinsX()+1); if (ov>0) cout << "WARNING: " << ov << " overflows in residuals." << endl; cout << dec << endl; cout << " Number of calls to FCN: " << minuit.fNfcn << endl; cout << "Minimum value found for FCN (Chi^2?): " << minuit.fAmin << endl; cout << " Fit-Probability(?): " << TMath::Prob(minuit.fAmin/*fOriginal.GetSize()*/, fOriginal.GetSize()-minuit.GetNumFreePars())*100 << "%" << endl; cout << " Chi^2/NDF: " << minuit.fAmin/(fOriginal.GetSize()-minuit.GetNumFreePars()) << endl; //cout << "Prob(?): " << TMath::Prob(fChisquare,ndf); // // Print all data sets for which the backward correction is // twice times worse than the residual gotten from the // bending correction itself // cout << endl; cout << "Checking backward correction (raw-->star):" << endl; for (int i=0; iFindObject("CanvGraphs"))) delete c1; if ((c1 = (TCanvas*)gROOT->FindObject("CanvResiduals"))) delete c1; c1=new TCanvas("CanvGraphs", "Graphs"); c1->Divide(2,3,0,0); c1->cd(1); gPad->SetBorderMode(0); TGraph *g=(TGraph*)gaz.DrawClone("A*"); g->SetBit(kCanDelete); g->GetHistogram()->SetXTitle("Az [\\circ]"); g->GetHistogram()->SetYTitle("\\Delta Az [\\circ]"); c1->cd(2); gPad->SetBorderMode(0); g=(TGraph*)gdaz.DrawClone("A*"); g->SetBit(kCanDelete); g->GetHistogram()->SetXTitle("Zd [\\circ]"); g->GetHistogram()->SetYTitle("\\Delta Az [\\circ]"); cout << "Mean dAz: " << g->GetMean(2) << " \xb1 " << g->GetRMS(2) << endl; c1->cd(3); gPad->SetBorderMode(0); g=(TGraph*)gdzd.DrawClone("A*"); g->SetBit(kCanDelete); g->GetHistogram()->SetXTitle("Az [\\circ]"); g->GetHistogram()->SetYTitle("\\Delta Zd [\\circ]"); cout << "Mean dZd: " << g->GetMean(2) << " \xb1 " << g->GetRMS(2) << endl; cout << endl; c1->cd(4); gPad->SetBorderMode(0); g=(TGraph*)gzd.DrawClone("A*"); g->SetBit(kCanDelete); g->GetHistogram()->SetXTitle("Zd [\\circ]"); g->GetHistogram()->SetYTitle("\\Delta Zd [\\circ]"); c1->cd(5); gPad->SetBorderMode(0); g=(TGraph*)graz.DrawClone("A*"); g->SetBit(kCanDelete); g->GetHistogram()->SetXTitle("Az [\\circ]"); g->GetHistogram()->SetYTitle("\\Delta [\\circ]"); c1->cd(6); gPad->SetBorderMode(0); g=(TGraph*)grzd.DrawClone("A*"); g->SetBit(kCanDelete); g->GetHistogram()->SetXTitle("Zd [\\circ]"); g->GetHistogram()->SetYTitle("\\Delta [\\circ]"); // // Print out the residual before and after correction in several // units // cout << fCoordinates.GetSize() << " data sets." << endl << endl; cout << "Total Spread of Residual:" << endl; cout << "-------------------------" << endl; cout << "before: " << hres1.GetMean() << " \xb1 " << hres1.GetRMS() << " deg " << endl; cout << "after: " << hres2.GetMean() << " \xb1 " << hres2.GetRMS() << " deg " << endl; cout << "before: " << (int)(hres1.GetMean()*60+.5) << " \xb1 " << (int)(hres1.GetRMS()*60+.5) << " arcmin" << endl; cout << "after: " << (int)(hres2.GetMean()*60+.5) << " \xb1 " << (int)(hres2.GetRMS()*60+.5) << " arcmin" << endl; cout << "before: " << (int)(hres1.GetMean()*60*60/23.4+.5) << " \xb1 " << (int)(hres1.GetRMS()*60*60/23.4+.5) << " pix" << endl; cout << "after: " << (int)(hres2.GetMean()*60*60/23.4+.5) << " \xb1 " << (int)(hres2.GetRMS()*60*60/23.4+.5) << " pix" << endl; cout << "after: " << (int)(hres2.GetMean()*16384/360+.5) << " \xb1 " << (int)(hres2.GetRMS()*16384/360+.5) << " SE" << endl; cout << endl; cout << "backw: " << (int)(hres3.GetMean()*60+.5) << " \xb1 " << (int)(hres3.GetRMS()*60+.5) << " arcmin" << endl; cout << endl; // ± gStyle->SetOptStat(1110); gStyle->SetStatFormat("6.2g"); c1=new TCanvas("CanvResiduals", "Residuals", 800, 800); c1->Divide(2, 2, 0, 0); c1->cd(2); gPad->SetBorderMode(0); hres1.SetLineColor(kRed); hres1.DrawCopy(); c1->cd(4); gPad->SetBorderMode(0); hres2.SetLineColor(kBlue); TH1 *h=hres2.DrawCopy(); TF1 f("mygaus", "(gaus)", 0, 1); f.SetLineColor(kMagenta/*6*/); f.SetLineWidth(1); f.SetParameter(0, h->GetBinContent(1)); f.FixParameter(1, 0); f.SetParameter(2, h->GetRMS()); h->Fit("mygaus", "QR"); hres3.SetLineColor(kGreen); hres3.SetLineStyle(kDashed); hres3.DrawCopy("same"); cout << "Gaus-Fit Sigma: " << f.GetParameter(2) << "\xb0" << endl; cout << "Fit-Probability: " << f.GetProb()*100 << "%" << endl; cout << " Chi^2/NDF: " << f.GetChisquare()/f.GetNDF() << endl; c1->cd(1); gPad->SetBorderMode(0); gPad->SetTheta(90); gPad->SetPhi(90); TH2F h2res1("Res2D1", " Dataset positions on the sky ", 32, 0, 360, 10, 0, 90); h2res1.SetBit(TH1::kNoStats); h2res1.DrawCopy("surf1pol"); gPad->Modified(); gPad->Update(); for (int i=0; icd(3); gPad->SetBorderMode(0); gPad->SetTheta(90); gPad->SetPhi(90); h2res1.SetTitle(" Arb. Residuals after correction (scaled) "); h2res1.DrawCopy("surf1pol"); gPad->Modified(); gPad->Update(); for (int i=0; i> set; // Read data from file [deg], it is stored in [rad] if (!fin) break; fOriginal.Add(new Set(set)); } cout << "Found " << fOriginal.GetSize()-size << " sets of coordinates "; cout << "(Total=" << fOriginal.GetSize() << ")" << endl; } Bool_t ProcessMessage(Long_t msg, Long_t mp1, Long_t mp2) { // cout << "Msg: " << hex << GET_MSG(msg) << endl; // cout << "SubMsg: " << hex << GET_SUBMSG(msg) << dec << endl; switch (GET_MSG(msg)) { case kC_COMMAND: switch (GET_SUBMSG(msg)) { case kCM_BUTTON: switch (mp1) { case kTbFit: Fit(); DisplayBending(); return kTRUE; case kTbLoad: fBending.Load("bending_magic.txt"); DisplayBending(); return kTRUE; case kTbSave: fBending.Save("bending_magic.txt"); return kTRUE; case kTbLoadStars: LoadStars(); DisplayData(); return kTRUE; case kTbReset: fBending.Clear(); DisplayBending(); return kTRUE; case kTbResetStars: fOriginal.Delete(); DisplayData(); return kTRUE; } return kTRUE; } return kTRUE; } return kTRUE; } void AddTextButton(TGCompositeFrame *f, TString txt, Int_t id=-1, TGLayoutHints *h=0) { TGButton *but = new TGTextButton(f, txt, id); but->Associate(this); f->AddFrame(but, h); fList->Add(but); } void AddCheckButton(TGCompositeFrame *f, TString txt, Int_t id=-1, TGLayoutHints *h=0) { TGButton *but = new TGCheckButton(f, txt, id); but->Associate(this); f->AddFrame(but, h); fList->Add(but); } TGLabel *AddLabel(TGCompositeFrame *f, TString txt, TGLayoutHints *h=0) { TGLabel *l = new TGLabel(f, txt/*, TGLabel::GetDefaultGC()(), fFont*/); f->AddFrame(l, h); fList->Add(l); fLabel.Add(l); return l; } void DisplayBending() { TArrayD par, err; fBending.GetParameters(par); fBending.GetError(err); TGLabel *l; for (int i=0; iSetText(Form("%.4f\xb0", par[i])); l = (TGLabel*)fLabel.At(MBending::GetNumPar()+i); l->SetText(Form("\xb1 %8.4f\xb0", err[i])); } } void DisplayData() { TGLabel *l = (TGLabel*)fLabel.At(3*MBending::GetNumPar()); l->SetText(Form("%d data sets loaded.", fOriginal.GetSize())); } public: ~MFit() { if (fFont) gVirtualX->DeleteFont(fFont); } MFit() : TGMainFrame(gClient->GetRoot(), 740, 370, kHorizontalFrame) { fCoordinates.SetOwner(); fOriginal.SetOwner(); fList = new MGList; fList->SetOwner(); fFont = gVirtualX->LoadQueryFont("7x13bold"); TGLayoutHints *hints0 = new TGLayoutHints(kLHintsExpandY, 7, 5, 5, 6); TGLayoutHints *hints1 = new TGLayoutHints(kLHintsExpandX|kLHintsExpandY, 5, 7, 5, 6); TGLayoutHints *hints2 = new TGLayoutHints(kLHintsCenterX|kLHintsCenterY, 5, 5, 5, 5); fList->Add(hints0); fList->Add(hints1); fList->Add(hints2); TGGroupFrame *grp1 = new TGGroupFrame(this, "Control", kVerticalFrame); AddFrame(grp1, hints0); fList->Add(grp1); TGGroupFrame *grp2 = new TGGroupFrame(this, "Parameters", kHorizontalFrame); AddFrame(grp2, hints1); fList->Add(grp2); TGLayoutHints *hints4 = new TGLayoutHints(kLHintsExpandX, 5, 5, 5); TGLayoutHints *hints5 = new TGLayoutHints(kLHintsExpandX, 5, 5, 15); AddTextButton(grp1, "Load Pointing Model", kTbLoad, hints5); AddTextButton(grp1, "Save Pointing Model", kTbSave, hints4); AddTextButton(grp1, "Fit Parameters", kTbFit, hints5); AddTextButton(grp1, "Reset Parameters", kTbReset, hints4); AddTextButton(grp1, "Load Stars", kTbLoadStars, hints5); AddTextButton(grp1, "Reset Stars", kTbResetStars, hints4); fList->Add(hints4); fList->Add(hints5); TGHorizontalFrame *comp = new TGHorizontalFrame(grp2, 1, 1); grp2->AddFrame(comp); fList->Add(comp); TGLayoutHints *hints3 = new TGLayoutHints(kLHintsLeft|kLHintsTop, 0, 20, 5, 0); fList->Add(hints3); TGLabel *l; TGVerticalFrame *vframe = new TGVerticalFrame(comp, 1, 1); comp->AddFrame(vframe, hints3); fList->Add(vframe); for (int i=0; iAddFrame(vframe, hints3); fList->Add(vframe); l = new TGLabel(vframe, "+000.0000"); l->SetTextJustify(kTextRight); fList->Add(l); fLabel.Add(l); TGButton *but = (TGButton*)fList->FindWidget(0); TGLayoutHints *h = new TGLayoutHints(kLHintsCenterY, 0, 0, but->GetHeight()-l->GetHeight()); fList->Add(h); vframe->AddFrame(l,h); for (int i=1; iSetTextJustify(kTextRight); vframe = new TGVerticalFrame(comp, 1, 1); comp->AddFrame(vframe, hints3); fList->Add(vframe); for (int i=0; iSetTextJustify(kTextRight); vframe = new TGVerticalFrame(comp, 1, 1); comp->AddFrame(vframe, hints3); fList->Add(vframe); for (int i=0; iAddFrame(l, hints5); fList->Add(l); fLabel.Add(l); ((TGCheckButton*)fList->FindWidget(0))->SetState(kButtonDown); ((TGCheckButton*)fList->FindWidget(1))->SetState(kButtonDown); ((TGCheckButton*)fList->FindWidget(3))->SetState(kButtonDown); ((TGCheckButton*)fList->FindWidget(4))->SetState(kButtonDown); ((TGCheckButton*)fList->FindWidget(5))->SetState(kButtonDown); ((TGCheckButton*)fList->FindWidget(6))->SetState(kButtonDown); ((TGCheckButton*)fList->FindWidget(7))->SetState(kButtonDown); SetWindowName("TPoint Fitting Window"); SetIconName("TPoint++"); Layout(); MapSubwindows(); MapWindow(); DisplayBending(); DisplayData(); } ClassDef(MFit, 0) }; ClassImp(MFit); void gui() { new MFit; }