Index: trunk/MagicSoft/Cosy/Changelog =================================================================== --- trunk/MagicSoft/Cosy/Changelog (revision 2406) +++ trunk/MagicSoft/Cosy/Changelog (revision 2407) @@ -1,3 +1,34 @@ -*-*- END -*-*- + 2003/10/20 - Thomas Bretz + + * base/File.cc: + - only close file if f!=NULL + + * base/MThread.cc: + - added debug output (depending on DEBUG macro) + + * candrv/vmodican.cc: + - do not call exit in case the module couldn't be opened + - changes some comments + + * catalog/StarCatalog.cc: + - do not exit in case the catalog was not found + - initialize fSao and fSrt with NULL + - only delete fSao and fSrt if != NULL + - only calculate stars if catalog was loaded successfully + + * main/MBending.[h,cc]: + - added file header + - added class decription + - added FLOP, TF and TX + - changed such that adding new stuff is easier + - changed order in correction + - replaced all switch statement by simple loops + + * tpoint/gui.C: + - full support for fitting a pointing correction + + + 2003/10/15 - Thomas Bretz (La Palma) Index: trunk/MagicSoft/Cosy/base/File.cc =================================================================== --- trunk/MagicSoft/Cosy/base/File.cc (revision 2406) +++ trunk/MagicSoft/Cosy/base/File.cc (revision 2407) @@ -12,5 +12,6 @@ File::~File() { - fclose(f); + if (f) + fclose(f); } Index: trunk/MagicSoft/Cosy/base/MThread.cc =================================================================== --- trunk/MagicSoft/Cosy/base/MThread.cc (revision 2406) +++ trunk/MagicSoft/Cosy/base/MThread.cc (revision 2407) @@ -7,4 +7,5 @@ #undef DEBUG +//#define DEBUG // ---------------------------------------------------------------------- @@ -91,5 +92,13 @@ fStop = false; +#ifdef DEBUG + cout << "MThread::RunThread" << endl; +#endif + void *rc = Thread(); + +#ifdef DEBUG + cout << "MThread::RunThread...done." << endl; +#endif fIsRunning = false; @@ -106,5 +115,7 @@ { MThread *thread = (MThread*)arg; - +#ifdef DEBUG + cout << "MThread::MapThread" << endl; +#endif return thread->RunThread(); } Index: trunk/MagicSoft/Cosy/candrv/vmodican.cc =================================================================== --- trunk/MagicSoft/Cosy/candrv/vmodican.cc (revision 2406) +++ trunk/MagicSoft/Cosy/candrv/vmodican.cc (revision 2407) @@ -762,6 +762,6 @@ // /* can_send - send message to standard host interface (mid priority) // * -// * This function sends a complete message to the standard host interface of -// * a VMOD-ICAN. +// * This function sends a complete message to the standard host interface +// * of a VMOD-ICAN. // * // * The message will be queued to the middle prioritized of the three @@ -951,5 +951,4 @@ cout << strerror(errno) << endl; return; -// exit(1); // open module } Index: trunk/MagicSoft/Cosy/catalog/StarCatalog.cc =================================================================== --- trunk/MagicSoft/Cosy/catalog/StarCatalog.cc (revision 2406) +++ trunk/MagicSoft/Cosy/catalog/StarCatalog.cc (revision 2407) @@ -14,5 +14,5 @@ ClassImp(StarCatalog); -StarCatalog::StarCatalog(MObservatory::LocationName_t key) : SlaStars(key), fEntries(0), fSinAngle(0), fCosAngle(1) +StarCatalog::StarCatalog(MObservatory::LocationName_t key) : SlaStars(key), fSao(NULL), fSrt(NULL), fEntries(0), fSinAngle(0), fCosAngle(1) { // p = pointer to MainFrame (not owner) @@ -22,7 +22,6 @@ // File idx("sao/sao-sort.idx", "r"); - if (!idx) - exit(0); + return; while (!idx.Eof()) @@ -52,6 +51,8 @@ StarCatalog::~StarCatalog() { - delete fSrt; - delete fSao; + if (fSrt) + delete fSrt; + if (fSao) + delete fSao; } @@ -647,4 +648,7 @@ // --------- search for stars in catalog ---------- // + if (fEntries==0) + return; + int count = 0; int deleted = 0; Index: trunk/MagicSoft/Cosy/main/MBending.cc =================================================================== --- trunk/MagicSoft/Cosy/main/MBending.cc (revision 2406) +++ trunk/MagicSoft/Cosy/main/MBending.cc (revision 2407) @@ -1,2 +1,60 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Thomas Bretz, 2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MBending +// +// Double_t fIe ; // [rad] Index Error in Elevation +// Double_t fIa ; // [rad] Index Error in Azimuth +// +// Double_t fFlop ; // [rad] Vertical Sag +// * do not use if not data: Zd<0 +// +// Double_t fNpae ; // [rad] Az-El Nonperpendicularity +// +// Double_t fCa ; // [rad] Left-Right Collimation Error +// +// Double_t fAn ; // [rad] Azimuth Axis Misalignment (N-S) +// Double_t fAw ; // [rad] Azimuth Axis Misalignment (E-W) +// +// Double_t fTf ; // [rad] Tube fluxture (sin) +// * same as ecec if no data: Zd<0 +// Double_t fTx ; // [rad] Tube fluxture (tan) +// * do not use with NPAE if no data: Zd<0 +// +// Double_t fNrx ; // [rad] Nasmyth rotator displacement, horizontan +// Double_t fNry ; // [rad] Nasmyth rotator displacement, vertical +// +// Double_t fCrx ; // [rad] Alt/Az Coude Displacement (N-S) +// Double_t fCry ; // [rad] Alt/Az Coude Displacement (E-W) +// +// Double_t fEces ; // [rad] Elevation Centering Error (sin) +// Double_t fAces ; // [rad] Azimuth Centering Error (sin) +// Double_t fEcec ; // [rad] Elevation Centering Error (cos) +// Double_t fAcec ; // [rad] Azimuth Centering Error (cos) +// +//////////////////////////////////////////////////////////////////////////// #include "MBending.h" @@ -10,4 +68,52 @@ ClassImp(MBending); +const Int_t MBending::fNumPar=19; + +void MBending::Init() +{ + fCoeff = new Double_t*[fNumPar]; + fName = new TString[fNumPar]; + fDescr = new TString[fNumPar]; + + fCoeff[ 0] = &fIa; fName[ 0] = "IA"; + fCoeff[ 1] = &fIe; fName[ 1] = "IE"; + fCoeff[ 2] = &fFlop; fName[ 2] = "FLOP"; + fCoeff[ 3] = &fNpae; fName[ 3] = "NPAE"; + fCoeff[ 4] = &fCa; fName[ 4] = "CA"; + fCoeff[ 5] = &fAn; fName[ 5] = "AN"; + fCoeff[ 6] = &fAw; fName[ 6] = "AW"; + fCoeff[ 7] = &fTf; fName[ 7] = "TF"; + fCoeff[ 8] = &fTx; fName[ 8] = "TX"; + fCoeff[ 9] = &fEces; fName[ 9] = "ECES"; + fCoeff[10] = &fAces; fName[10] = "ACES"; + fCoeff[11] = &fEcec; fName[11] = "ECEC"; + fCoeff[12] = &fAcec; fName[12] = "ACEC"; + fCoeff[13] = &fNrx; fName[13] = "NRX"; + fCoeff[14] = &fNry; fName[14] = "NRY"; + fCoeff[15] = &fCrx; fName[15] = "CRX"; + fCoeff[16] = &fCry; fName[16] = "CRY"; + fCoeff[17] = &fMagic1; fName[17] = "MAGIC1"; + fCoeff[18] = &fMagic2; fName[18] = "MAGIC2"; + + fDescr[ 0] = "Index Error Azimuth"; + fDescr[ 1] = "Index Error Zenith Distance"; + fDescr[ 2] = "Vertical Sag"; + fDescr[ 3] = "Az-El Nonperpendicularity"; + fDescr[ 4] = "Left-Right Collimation Error"; + fDescr[ 5] = "Azimuth Axis Misalignment (N-S)"; + fDescr[ 6] = "Azimuth Axis Misalignment (E-W)"; + fDescr[ 7] = "Tube fluxture (sin)"; + fDescr[ 8] = "Tube fluxture (tan)"; + fDescr[ 9] = "Elevation Centering Error (sin)"; + fDescr[10] = "Azimuth Centering Error (sin)"; + fDescr[11] = "Elevation Centering Error (cos)"; + fDescr[12] = "Azimuth Centering Error (cos)"; + fDescr[13] = "Nasmyth rotator displacement (horizontal)"; + fDescr[14] = "Nasmyth rotator displacement (vertical)"; + fDescr[15] = "Alt/Az Coude Displacement (N-S)"; + fDescr[16] = "Alt/Az Coude Displacement (E-W)"; + fDescr[17] = "n/a"; + fDescr[18] = "n/a"; +} void MBending::Reset() { @@ -62,6 +168,6 @@ cout << endl; - cout << " & = Name Value Sigma" << endl; - cout << "------------------------------------------" << endl; + cout << " & = Name Value Sigma" << endl; + cout << "--------------------------------------------------" << endl; while (fin) @@ -90,26 +196,40 @@ fin >> val; - cout << str << "\t" << val << "° \t"; + cout << str << "\t" << setw(11) << val << "° \t"; val *= kDeg2Rad; - if (str=="IA") fIa = val; - if (str=="IE") fIe = val; - if (str=="CA") fCa = val; - if (str=="AN") fAn = val; - if (str=="AW") fAw = val; - if (str=="NRX") fNrx = val; - if (str=="NRY") fNry = val; - if (str=="CRX") fCrx = val; - if (str=="CRY") fCry = val; - if (str=="NPAE") fNpae = val; - if (str=="ECES") fEces = val; - if (str=="ACES") fAces = val; - if (str=="ECEC") fEcec = val; - if (str=="ACEC") fAcec = val; - //if (str=="MAGIC1") fMagic1 = val; - //if (str=="MAGIC2") fMagic2 = val; + Double_t *dest=NULL; + + if (str=="IA") dest = &fIa; + if (str=="IE") dest = &fIe; + if (str=="FLOP") dest = &fFlop; + if (str=="NPAE") dest = &fNpae; + if (str=="CA") dest = &fCa; + if (str=="AN") dest = &fAn; + if (str=="AW") dest = &fAw; + if (str=="TF") dest = &fTf; + if (str=="TX") dest = &fTx; + if (str=="NRX") dest = &fNrx; + if (str=="NRY") dest = &fNry; + if (str=="CRX") dest = &fCrx; + if (str=="CRY") dest = &fCry; + if (str=="ECES") dest = &fEces; + if (str=="ACES") dest = &fAces; + if (str=="ECEC") dest = &fEcec; + if (str=="ACEC") dest = &fAcec; + + if (dest) + *dest = val; fin >> val; - cout << val*kDeg2Rad << endl; + cout << setw(9) << val << "°" << endl; + + // Find corresponding error + for (int i=0; i +#endif class TMinuit; @@ -9,11 +13,16 @@ { private: + static const Int_t fNumPar; + Double_t fIe ; // [rad] Index Error in Elevation Double_t fIa ; // [rad] Index Error in Azimuth + Double_t fFlop ; // [rad] Vertical Sag + Double_t fNpae ; // [rad] Az-El Nonperpendicularity Double_t fCa ; // [rad] Left-Right Collimation Error Double_t fAn ; // [rad] Azimuth Axis Misalignment (N-S) Double_t fAw ; // [rad] Azimuth Axis Misalignment (E-W) - Double_t fNpae ; // [rad] Az-El Nonperpendicularity - Double_t fNrx ; // [rad] Nasmyth rotator displacement, horizontan + Double_t fTf ; // [rad] Tube fluxture (sin) + Double_t fTx ; // [rad] Tube fluxture (tan) + Double_t fNrx ; // [rad] Nasmyth rotator displacement, horizontal Double_t fNry ; // [rad] Nasmyth rotator displacement, vertical Double_t fCrx ; // [rad] Alt/Az Coude Displacement (N-S) @@ -26,28 +35,27 @@ Double_t fMagic2; // [rad] Magic Term (what is it?) + Double_t **fCoeff; //! + TString *fName; //! + TString *fDescr; //! + + TArrayD fError; + + void Init(); void Clear() { - fIe =0 ; // Index Error in Elevation - fIa =0 ; // Index Error in Azimuth - fCa =0 ; // Left-Right Collimation Error - fAn =0 ; // Azimuth Axis Misalignment (N-S) - fAw =0 ; // Azimuth Axis Misalignment (E-W) - fNpae=0 ; // Az-El Nonperpendicularity - fNrx =0 ; // Nasmyth rotator displacement, horizontan - fNry =0 ; // Nasmyth rotator displacement, vertical - fCrx =0 ; // Alt/Az Coude Displacement (N-S) - fCry =0 ; // Alt/Az Coude Displacement (E-W) - fEces=0 ; // Elevation Centering Error (sin) - fAces=0 ; // Azimuth Centering Error (sin) - fEcec=0 ; // Elevation Centering Error (cos) - fAcec=0 ; // Azimuth Centering Error (cos) - fMagic1=0; // Azimuth Centering Error (cos) - fMagic2=0; // Azimuth Centering Error (cos) + for (int i=0; i #include +#include #include @@ -5,6 +7,121 @@ #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 @@ -12,22 +129,553 @@ private: enum { - kCbIA, - kCbIE, - kCbNPAE, - kCbCA, - kCbAN, - kCbAW, - kCbCRX, - kCbCRY, - kCbNRX, - kCbNRY + kTbFit = 19, //MBending::GetNumPar(), // FIXME!!! + kTbLoad, + kTbSave, + kTbLoadStars, + kTbReset, + kTbResetStars }; - MGList *list; + 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) << endl; + // cout << "Msg: " << hex << GET_MSG(msg) << endl; + // cout << "SubMsg: " << hex << GET_SUBMSG(msg) << dec << endl; switch (GET_MSG(msg)) { @@ -35,9 +683,31 @@ switch (GET_SUBMSG(msg)) { - case kCM_CHECKBUTTON: - /* - TGCheckButton *but = (TGCheckButton*)list->FindWidget(mp1); - cout << hex << but->GetName() << " " << mp2 << endl; - */ + 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; } @@ -52,5 +722,5 @@ but->Associate(this); f->AddFrame(but, h); - list->Add(but); + fList->Add(but); } @@ -61,53 +731,91 @@ but->Associate(this); f->AddFrame(but, h); - list->Add(but); - } - - void AddLabel(TGCompositeFrame *f, TString txt, TGLayoutHints *h=0, TList *label=0) - { - TGLabel *l = new TGLabel(f, txt); + 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); - list->Add(l); - if (label) - label->Add(l); + 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() : TGMainFrame(gClient->GetRoot(), 550, 350, kHorizontalFrame) - { - list = new MGList; - list->SetOwner(); + ~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); - list->Add(hints0); - list->Add(hints1); - list->Add(hints2); + fList->Add(hints0); + fList->Add(hints1); + fList->Add(hints2); TGGroupFrame *grp1 = new TGGroupFrame(this, "Control", kVerticalFrame); AddFrame(grp1, hints0); - list->Add(grp1); + fList->Add(grp1); TGGroupFrame *grp2 = new TGGroupFrame(this, "Parameters", kHorizontalFrame); AddFrame(grp2, hints1); - list->Add(grp2); - - - - TGLayoutHints *hints4 = new TGLayoutHints(kLHintsExpandX, 5, 5, 10); - AddTextButton(grp1, "Load Bending Model", -1, hints4); - AddTextButton(grp1, "Save Bending Model", -1, hints4); - AddTextButton(grp1, "Load Observations", -1, hints4); - AddTextButton(grp1, "Fit Parameters", -1, hints4); - AddTextButton(grp1, "Reset Parameters", -1, hints4); + 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); - list->Add(comp); + fList->Add(comp); TGLayoutHints *hints3 = new TGLayoutHints(kLHintsLeft|kLHintsTop, 0, 20, 5, 0); - list->Add(hints3); + fList->Add(hints3); TGLabel *l; @@ -115,58 +823,58 @@ TGVerticalFrame *vframe = new TGVerticalFrame(comp, 1, 1); comp->AddFrame(vframe, hints3); - list->Add(vframe); - - AddCheckButton(vframe, "IA", kCbIA); - AddCheckButton(vframe, "IE", kCbIE); - AddCheckButton(vframe, "NPAE", kCbNPAE); - AddCheckButton(vframe, "CA", kCbCA); - AddCheckButton(vframe, "AN", kCbAN); - AddCheckButton(vframe, "AW", kCbAW); - AddCheckButton(vframe, "CRX", kCbCRX); - AddCheckButton(vframe, "CRY", kCbCRY); - AddCheckButton(vframe, "NRX", kCbNRX); - AddCheckButton(vframe, "NRY", kCbNRY); + fList->Add(vframe); + + for (int i=0; iAddFrame(vframe, hints3); - list->Add(vframe); - - - TList label; - - l = new TGLabel(vframe, "0"); - list->Add(l); - label.Add(l); - - TGButton *but = (TGButton*)list->FindWidget(kCbIA); + 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()); - list->Add(h); + fList->Add(h); vframe->AddFrame(l,h); - AddLabel(vframe, "0", h, &label); - AddLabel(vframe, "0", h, &label); - AddLabel(vframe, "0", h, &label); - AddLabel(vframe, "0", h, &label); - AddLabel(vframe, "0", h, &label); - AddLabel(vframe, "0", h, &label); - AddLabel(vframe, "0", h, &label); - AddLabel(vframe, "0", h, &label); - AddLabel(vframe, "0", h, &label); + for (int i=1; iSetTextJustify(kTextRight); vframe = new TGVerticalFrame(comp, 1, 1); comp->AddFrame(vframe, hints3); - list->Add(vframe); - - AddLabel(vframe, "Offset Azimuth", h, &label); - AddLabel(vframe, "Offset Zenith Distance", h, &label); - AddLabel(vframe, "0", h, &label); - AddLabel(vframe, "0", h, &label); - AddLabel(vframe, "0", h, &label); - AddLabel(vframe, "0", h, &label); - AddLabel(vframe, "0", h, &label); - AddLabel(vframe, "0", h, &label); - AddLabel(vframe, "0", h, &label); + 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(); @@ -175,16 +883,11 @@ MapWindow(); - - - - l = (TGLabel*)label.At(0); - l->SetText("221.4567"); - - l = (TGLabel*)label.At(1); - l->SetText("-1.65423"); - - Layout(); - } + DisplayBending(); + DisplayData(); + } + ClassDef(MFit, 0) }; + +ClassImp(MFit); void gui() @@ -192,3 +895,2 @@ new MFit; } -