Index: /trunk/MagicSoft/Cosy/Changelog =================================================================== --- /trunk/MagicSoft/Cosy/Changelog (revision 1809) +++ /trunk/MagicSoft/Cosy/Changelog (revision 1810) @@ -1,3 +1,48 @@ -*-*- END -*-*- + 2003/03/11 - Daniela Dorner, Thomas Bretz: + + * base/MStar.h: + - added Compare + + * base/MStarList.cc: + - some small bugfixes + + * base/MStarList.h: + - added Sort + - added Expand + + * base/timer.[h,cc]: + - Added GetTimeval + + * gui/MGAccuracy.cc: + - Exchanged zd/az in calculation of Residual!!! + + * gui/MGPngReader.[h,cc]: + - set default lim mag to 7.0 + - added new ouput for the pointing position fPZdAz + - added/fixed TrackingError/CalcTrackingError + - changed Filter2 to CaosFilter + - reordered starguider stuff in Execute + - changed color of circles + + * main/MBending.[h,cc]: + - removed MAGIC1 and MAGIC2 + - removed '-' from writing + - fixed some bugs in the enumerations of the coefficients + - added some formating option for output + + * tpoint/tpointfit.C: + - removed usage of MyAdjust + - fixed the Calculation of the residuals + - fixed reading + - added some correction in case of an overflow (360deg/0deg) + - fixed drawing + - added second Migrad turn... + - changed the screen and graphical output + + * videodev/CaosFilter.[h,cc]: + - changed RemoveTwins to accept a radius + + 2003/03/02 - Daniela Dorner, Thomas Bretz (LaPalma): Index: /trunk/MagicSoft/Cosy/base/MStar.h =================================================================== --- /trunk/MagicSoft/Cosy/base/MStar.h (revision 1809) +++ /trunk/MagicSoft/Cosy/base/MStar.h (revision 1810) @@ -21,4 +21,20 @@ void Set(Double_t mx, Double_t my) { fX=mx; fY=my; } + Int_t Compare(const TObject *obj) const + { + const MStar *const s = (MStar*)obj; + + if (fMagfMag) + return -1; + + if (fMag>s->fMag) + return 1; + + return 0; + } + + Bool_t IsSortable() const { return kTRUE; } + + ClassDef(MStar, 1) }; Index: /trunk/MagicSoft/Cosy/base/MStarList.cc =================================================================== --- /trunk/MagicSoft/Cosy/base/MStarList.cc (revision 1809) +++ /trunk/MagicSoft/Cosy/base/MStarList.cc (revision 1810) @@ -1,3 +1,5 @@ #include "MStarList.h" + +#include void MStarList::RemoveTwins(Double_t radius) @@ -18,6 +20,7 @@ return; + fStars.RemoveAt(idx); + MStarListIter Next(this, *first, radius); - Delete(idx); MStar *pos; @@ -29,5 +32,5 @@ mx += pos->GetX(); my += pos->GetY(); - Delete(pos); + fStars.Remove(pos); cnt++; } Index: /trunk/MagicSoft/Cosy/base/MStarList.h =================================================================== --- /trunk/MagicSoft/Cosy/base/MStarList.h (revision 1809) +++ /trunk/MagicSoft/Cosy/base/MStarList.h (revision 1810) @@ -38,4 +38,8 @@ void RemoveTwins(Double_t radius); + + void Sort() { fStars.Sort(); } + + void Expand(int n) { fStars.Expand(n); } }; Index: /trunk/MagicSoft/Cosy/base/timer.cc =================================================================== --- /trunk/MagicSoft/Cosy/base/timer.cc (revision 1809) +++ /trunk/MagicSoft/Cosy/base/timer.cc (revision 1810) @@ -26,4 +26,11 @@ // fDiv = fmod((fMs+fSecs)/(60*60*24), 1.0); } + +void Timer::GetTimeval(struct timeval *tv) const +{ + tv->tv_sec = fSecs; + tv->tv_usec = fMs; +} + /* void Timer::Set(const long mjd) Index: /trunk/MagicSoft/Cosy/base/timer.h =================================================================== --- /trunk/MagicSoft/Cosy/base/timer.h (revision 1809) +++ /trunk/MagicSoft/Cosy/base/timer.h (revision 1810) @@ -33,4 +33,6 @@ void SetTimer(const struct timeval *tv); + void GetTimeval(struct timeval *tv) const; + int GetSecs() { return fSecs; } Index: /trunk/MagicSoft/Cosy/gui/MGAccuracy.cc =================================================================== --- /trunk/MagicSoft/Cosy/gui/MGAccuracy.cc (revision 1809) +++ /trunk/MagicSoft/Cosy/gui/MGAccuracy.cc (revision 1810) @@ -195,8 +195,10 @@ aaz *= d2r; + const double el = TMath::Pi()/2-pzd; + const double dphi2 = aaz/2.; const double cos2 = cos(dphi2)*cos(dphi2); const double sin2 = sin(dphi2)*sin(dphi2); - const double d = cos(azd)*cos2 - cos(2*pzd+azd)*sin2; + const double d = cos(azd)*cos2 - cos(2*el+azd)*sin2; double dist = acos(d); @@ -250,5 +252,5 @@ UpdateCross(x, y); - UpdateText(pos.Zd(), acc.Az(), acc.Zd()); + UpdateText(pos.Zd(), acc.Zd(), acc.Az()); SetModified(); Index: /trunk/MagicSoft/Cosy/gui/MGPngReader.cc =================================================================== --- /trunk/MagicSoft/Cosy/gui/MGPngReader.cc (revision 1809) +++ /trunk/MagicSoft/Cosy/gui/MGPngReader.cc (revision 1810) @@ -3,8 +3,10 @@ #include // ifstream #include // cout +#include // cout #include #include +#include #include #include // TGHorizontal3DLine @@ -13,4 +15,5 @@ #include "MGImage.h" #include "MGCoordinates.h" +#include "MGDispCoordinates.h" #include "coord.h" @@ -255,5 +258,5 @@ fList->Add(fLimMag); - fSao->SetLimitMag(8.0); + fSao->SetLimitMag(7.0); fInterpol = new MGPopupMenu(p); @@ -314,9 +317,32 @@ fList->Add(fCRaDec); - fCZdAz = new MGCoordinates(this, kETypeZdAz); + fCZdAz = new MGCoordinates(this, kETypeZdAz, kFALSE); fCZdAz->Move(240+12+10, fMenu->GetDefaultHeight()+584); AddFrame(fCZdAz); fList->Add(fCZdAz); + fPZdAz = new MGCoordinates(this, kETypeZdAz, kFALSE); + fPZdAz->Move(240+12+10, fMenu->GetDefaultHeight()+630); + AddFrame(fPZdAz); + fList->Add(fPZdAz); + + TGLabel *l = new TGLabel(this, "Arb.-Sky Pos"); + l->SetTextJustify(kTextLeft); + l->Move(480+32, fMenu->GetDefaultHeight()+590); + AddFrame(l); + fList->Add(l); + + l = new TGLabel(this, "arcsec/pix"); + l->SetTextJustify(kTextLeft); + l->Move(605, fMenu->GetDefaultHeight()+619); + AddFrame(l); + fList->Add(l); + + l = new TGLabel(this, "Pointing Pos"); + l->SetTextJustify(kTextLeft); + l->Move(480+32, fMenu->GetDefaultHeight()+655); + AddFrame(l); + fList->Add(l); + const Double_t pixsize = 23.4; @@ -328,5 +354,5 @@ fPixSize = new TGTextEntry(this, txt, IDM_kPixSize); fPixSize->SetAlignment(kTextCenterX); - fPixSize->Move(600, fMenu->GetDefaultHeight()+584); + fPixSize->Move(547, fMenu->GetDefaultHeight()+617); AddFrame(fPixSize); fList->Add(fPixSize); @@ -357,8 +383,14 @@ MapSubwindows(); MapWindow(); + + //------------------------------------------------------------ + // XY xy(3.819444, 24.05333); + // fCRaDec->SetCoordinates(xy); + // fRaDec->Set(xy.X()*360/24, xy.Y()); + //------------------------------------------------------------ } MGPngReader::MGPngReader(MObservatory::LocationName_t obs) -: TGMainFrame(gClient->GetRoot(), 768, 700), fFile(NULL), fDx((768-kZOOM)/2), fDy((512-kZOOM)/2) +: TGMainFrame(gClient->GetRoot(), 768, 740), fFile(NULL), fDx((768-kZOOM)/2), fDy((512-kZOOM)/2) { fSao = new StarCatalog(obs); @@ -811,6 +843,71 @@ } -void MGPngReader::CalcTrackingError(MStarList &spots, MStarList &stars) -{ +ZdAz MGPngReader::TrackingError(TArrayF &x, TArrayF &y, TArrayF &mag) const +{ + // + // Viewable area (FIXME: AZ) + // + TH2F h("Hist", "dX/dY", 77, -768/2-.5, 768/2+.5, 58, -576/2-.5, 576/2+.5); // 3 + + /* + TH1F hmag("HistMag", "Mag", 19, 0, 100); + for (int i=0; i48-15 && mag[i]<48+15) + h.Fill(x[i], y[i]); + } + + // + // Serach for an excess in the histogram + // + Int_t mx, my, dummy; + h.GetMaximumBin(mx, my, dummy); + + const double xmax = h.GetXaxis()->GetBinCenter(mx); + const double dx = h.GetXaxis()->GetBinWidth(mx); + + const double ymax = h.GetYaxis()->GetBinCenter(my); + const double dy = h.GetYaxis()->GetBinWidth(my); + + cout << setprecision(3); + cout << "Cut-XY: " << xmax << " +- " << dx << " / " << ymax << " +- " << dy << endl; + + TGraph g; + for (int i=0; ixmax-dx && x[i]ymax-dy && y[i]48-15 && mag[i]<48+15)) + continue; + + g.SetPoint(g.GetN(), x[i], y[i]); + } + + cout << "Offset-XY: " << g.GetMean(1) << " +- " << g.GetRMS(1) << " / "; + cout << g.GetMean(2) << " +- " << g.GetRMS(2) << endl; + + AltAz pos0 = fSao->CalcAltAzFromPix(768/2, 576/2)*kRad2Deg; + AltAz pos1 = fSao->CalcAltAzFromPix(768/2+g.GetMean(1), 576/2+g.GetMean(2))*kRad2Deg; + + pos1 -= pos0; + + ofstream fout("tracking_error.txt"); + fout << setprecision(10) << fSao->GetMjd()-52000 << " " << -pos1.Alt() << " " << pos1.Az() << endl; + fout.close(); + + return ZdAz(-pos1.Alt(), pos1.Az()); +} + +void MGPngReader::CalcTrackingError(Leds &leds, MStarList &stars) +{ + const Int_t max = leds.GetEntries(); + if (stars.GetRealEntries() < 3) { @@ -819,5 +916,5 @@ } - if (spots.GetRealEntries() < 1) + if (max < 1) { cout << "Sorry, less than 1 detected spot in FOV!" << endl; @@ -825,82 +922,59 @@ } - Int_t idx = 0; - - MStarList sortedspots; + stars.Sort(); // Sort by magnitude + + TString str = "data/tracking_"; + str += fSao->GetMjd()-52000; + str += ".txt"; + + ofstream fout(str); + + TArrayF x, y, mag; + + Int_t num = 0; + + // FIXME: Is predifined value 3 a good idea? MStar *star; - MStar *spot; MStarListIter NextStar(&stars); - MStarListIter NextSpot(&spots); - - while ((spot=NextSpot())) - { - AltAz aa = fSao->CalcAltAzFromPix(spot->GetX(), spot->GetY()); - spot->Set(aa.Az(), aa.Alt()); - } - - while ((star=NextStar())) - { - AltAz aa = fSao->CalcAltAzFromPix(star->GetX(), star->GetY()); - star->Set(aa.Az(), aa.Alt()); - - const double aaz = star->GetX(); - const double dphi2 = aaz/2.; - const double cos2 = cos(dphi2)*cos(dphi2); - const double sin2 = sin(dphi2)*sin(dphi2); - - Double_t min = 800; - - NextSpot.Reset(); - while ((spot=NextSpot())) + while ((star=NextStar()) && num++GetY(); - const double azd = TMath::Pi()/2-star->GetY(); - - const double d = cos(azd)*cos2 - cos(2*pzd+azd)*sin2; - - const Double_t dist = acos(d); - - if (dist>=min) - continue; - - min = dist; - sortedspots.AddAt(idx, spot->GetX(), spot->GetY(), spot->GetMag()); + const XY dpos(spot->GetX()-star->GetX(), spot->GetY()-star->GetY()); + + const Int_t idx = x.GetSize(); + + x.Set(idx+1); + y.Set(idx+1); + mag.Set(idx+1); + + x.AddAt(dpos.X(), idx); + y.AddAt(dpos.Y(), idx); + mag.AddAt(spot->GetMag()/star->GetMag(), idx); + + if (fout) + fout << x[idx] << " " << y[idx] << " " << mag[idx] << endl; } - if (min>768) - { - cout << "ERROR!!!!!!!!" << endl; - return; - } - idx++; - } - - // - // Now we have in sortedspots the entries with the shortest distances - // to the corresponding ones in stars. - // Now calculate the tracking error. - // - NextStar.Reset(); - MStarListIter NextSpot2(&sortedspots); - - Double_t meanx=0; - Double_t meany=0; - - while ((star=NextStar())) - { - spot = NextSpot2(); - - meanx += star->GetX() - spot->GetX(); - meany += star->GetY() - spot->GetY(); - } - - meanx /= idx; - meany /= idx; - - cout << "Tracking Error: dAlt=" << meany*180/TMath::Pi(); - cout << "° dAz=" << meanx*180/TMath::Pi() << "° (calculated"; - cout << " with " << idx << " stars/spots)" << endl; -} - + } + + ZdAz d = TrackingError(x, y, mag); + + // + // Calculated offsets + // + + // round= floor(x+.5) + cout << "Offset-ZdAz: " << d.Zd()*60 << "' / " << d.Az()*60 << "'" << endl; + cout << "Offset-ZdAz: " << d.Zd()/360*16384 << " / " << d.Az()/360*16384 << " (SE) " << endl; + + // + // Current Pointing position + // + ZdAz cpos = fSao->GetZdAz()-d; + fPZdAz->SetCoordinates(cpos); +} @@ -945,9 +1019,12 @@ MStarList spots; + + /* if (fDisplay->IsEntryChecked(IDM_kStarguider)) Filter2::Execute(spots, c); else - if (fDisplay->IsEntryChecked(IDM_kFilter)) - Filter::Execute(c); + */ + if (fDisplay->IsEntryChecked(IDM_kFilter)) + Filter::Execute(c); if (fDisplay->IsEntryChecked(IDM_kCaosFilter)) @@ -1031,25 +1108,36 @@ if (fDisplay->IsEntryChecked(IDM_kCatalog)) { + Timer time(tm); + + GetCoordinates(); + + MStarList stars; + fSao->GetStars(stars, time.GetMjd(), *fRaDec); + + if (fDisplay->IsEntryChecked(IDM_kStarguider)) + { + Leds leds; + CaosFilter::Execute(c, leds, 1); + + cout << "Found: " << leds.GetEntries() << " leds" << endl; + + CaosFilter::RemoveTwins(leds, 3); + leds.Compress(); + + cout << "Rest: " << leds.GetEntries() << " leds" << endl; + + CalcTrackingError(leds, stars); + } + + byte cimg[768*576]; + fSao->GetImg(c, cimg, stars); + DrawCircle(c, 0.5); DrawCircle(c, 1.0); DrawCircle(c, 1.5); - byte cimg[768*576]; - - GetCoordinates(); - - Timer time(tm); - - MStarList stars; - fSao->GetStars(stars, time.GetMjd(), *fRaDec); - fSao->GetImg(c, cimg, stars); - //fSao->GetImg(c, cimg, time.CalcMjd(), *fRaDec); + fCZdAz->SetCoordinates(fSao->GetZdAz()); fImage->DrawColImg(c, cimg); - - fCZdAz->SetCoordinates(fSao->GetZdAz()); - - if (fDisplay->IsEntryChecked(IDM_kStarguider)) - CalcTrackingError(spots, stars); } else @@ -1067,8 +1155,8 @@ { const int dy = (int)sqrt(rpix*rpix-dx*dx); - img[cx+dx + (cy-dy)*768] = 0xd0; - img[cx+dx + (cy+dy)*768] = 0xd0; - img[cx-dy + (cy+dx)*768] = 0xd0; - img[cx+dy + (cy+dx)*768] = 0xd0; + img[cx+dx + (cy-dy)*768] = 0x40; + img[cx+dx + (cy+dy)*768] = 0x40; + img[cx-dy + (cy+dx)*768] = 0x40; + img[cx+dy + (cy+dx)*768] = 0x40; } } Index: /trunk/MagicSoft/Cosy/gui/MGPngReader.h =================================================================== --- /trunk/MagicSoft/Cosy/gui/MGPngReader.h (revision 1809) +++ /trunk/MagicSoft/Cosy/gui/MGPngReader.h (revision 1810) @@ -13,11 +13,10 @@ #include "MGImage.h" -#include -#include -#include -#include +#include "coord.h" -class AltAz; -class RaDec; +class TArrayF; +class TH1F; +class TH2F; +class TGraph; class TTimer; @@ -87,4 +86,6 @@ MGCoordinates *fCZdAz; + MGCoordinates *fPZdAz; + TGTextEntry *fPixSize; @@ -104,5 +105,6 @@ void Toggle(MGPopupMenu *p, UInt_t id); void GetCoordinates(); - void CalcTrackingError(MStarList &, MStarList &); + void CalcTrackingError(Leds &, MStarList &); + ZdAz TrackingError(TArrayF &alt, TArrayF &az, TArrayF &mag) const; void InitHists(); Index: /trunk/MagicSoft/Cosy/main/MBending.cc =================================================================== --- /trunk/MagicSoft/Cosy/main/MBending.cc (revision 1809) +++ /trunk/MagicSoft/Cosy/main/MBending.cc (revision 1810) @@ -107,6 +107,6 @@ if (str=="ECEC") fEcec = val; if (str=="ACEC") fAcec = val; - if (str=="MAGIC1") fMagic1 = val; - if (str=="MAGIC2") fMagic2 = val; + //if (str=="MAGIC1") fMagic1 = val; + //if (str=="MAGIC2") fMagic2 = val; fin >> val; @@ -143,20 +143,20 @@ fout << "S 00 000000 000000 0000000" << endl; fout << setprecision(8); - fout << " IA " << -kRad2Deg*fIa << " -1" << endl; - fout << " IE " << -kRad2Deg*fIe << " -1" << endl; - fout << " CA " << -kRad2Deg*fNpae << " -1" << endl; - fout << " NPAE " << -kRad2Deg*fCa << " -1" << endl; - fout << " AN " << -kRad2Deg*fAn << " -1" << endl; - fout << " AW " << -kRad2Deg*fAw << " -1" << endl; - fout << " NRX " << -kRad2Deg*fNrx << " -1" << endl; - fout << " NRY " << -kRad2Deg*fNry << " -1" << endl; - fout << " CRX " << -kRad2Deg*fCrx << " -1" << endl; - fout << " CRY " << -kRad2Deg*fCry << " -1" << endl; - fout << " ECES " << -kRad2Deg*fEces << " -1" << endl; - fout << " ACES " << -kRad2Deg*fAces << " -1" << endl; - fout << " ECEC " << -kRad2Deg*fEcec << " -1" << endl; - fout << " ACEC " << -kRad2Deg*fAcec << " -1" << endl; - fout << " MAGIC1 " << -kRad2Deg*fMagic1 << " -1" << endl; - fout << " MAGIC2 " << -kRad2Deg*fMagic2 << " -1" << endl; + fout << " IA " << kRad2Deg*fIa << " -1" << endl; + fout << " IE " << kRad2Deg*fIe << " -1" << endl; + fout << " CA " << kRad2Deg*fNpae << " -1" << endl; + fout << " NPAE " << kRad2Deg*fCa << " -1" << endl; + fout << " AN " << kRad2Deg*fAn << " -1" << endl; + fout << " AW " << kRad2Deg*fAw << " -1" << endl; + fout << " NRX " << kRad2Deg*fNrx << " -1" << endl; + fout << " NRY " << kRad2Deg*fNry << " -1" << endl; + fout << " CRX " << kRad2Deg*fCrx << " -1" << endl; + fout << " CRY " << kRad2Deg*fCry << " -1" << endl; + fout << " ECES " << kRad2Deg*fEces << " -1" << endl; + fout << " ACES " << kRad2Deg*fAces << " -1" << endl; + fout << " ECEC " << kRad2Deg*fEcec << " -1" << endl; + fout << " ACEC " << kRad2Deg*fAcec << " -1" << endl; + //fout << " MAGIC1 " << -kRad2Deg*fMagic1 << " -1" << endl; + //fout << " MAGIC2 " << -kRad2Deg*fMagic2 << " -1" << endl; fout << "END" << endl; } @@ -173,4 +173,7 @@ p += CEC; + // const AltAz MAGIC(0, -fMagic1*tan(p.Alt())-fMagic2); + // p += MAGIC; + const AltAz CRX(-fCrx*sin(p.Az()-p.Alt()), fCrx*cos(p.Az()-p.Alt())/cos(p.Alt())); const AltAz CRY(-fCry*cos(p.Az()-p.Alt()), -fCry*sin(p.Az()-p.Alt())/cos(p.Alt())); @@ -183,6 +186,6 @@ p += NRY; - const AltAz MAGIC(-fMagic1*sin(p.Az()-fMagic2), 0); - p += MAGIC; + // const AltAz MAGIC(-fMagic1*sin(p.Az()-fMagic2), 0); + // p += MAGIC; const AltAz AW( fAw*sin(p.Az()), -fAw*cos(p.Az())*tan(p.Alt())); @@ -191,7 +194,6 @@ p += AN; -// const AltAz MAGIC(-fMagic1*sin(p.Az()-fMagic2), 0); -// p += MAGIC; - + // const AltAz MAGIC(-fMagic1*sin(p.Az()-fMagic2), 0); + // p += MAGIC; const AltAz CA(0, -fCa/cos(p.Alt())); @@ -227,6 +229,6 @@ p -= AW; - const AltAz MAGIC(-fMagic1*sin(p.Az()-fMagic2), 0); - p -= MAGIC; + // const AltAz MAGIC(-fMagic1*sin(p.Az()-fMagic2), 0); + // p -= MAGIC; const AltAz NRY(fNry*cos(p.Alt()), -fNry*tan(p.Alt())); @@ -285,23 +287,23 @@ fEces =par[10]/kRad2Deg; // Elevation Centering Error (sin) case 10: - fCry =par[9]/kRad2Deg; // Alt/Az Coude Displacement (E-W) + fCry =par[9] /kRad2Deg; // Alt/Az Coude Displacement (E-W) case 9: - fCrx =par[8]/kRad2Deg; // Alt/Az Coude Displacement (N-S) + fCrx =par[8] /kRad2Deg; // Alt/Az Coude Displacement (N-S) case 8: - fNry =par[7]/kRad2Deg; // Nasmyth rotator displacement, vertical + fNry =par[7] /kRad2Deg; // Nasmyth rotator displacement, vertical case 7: - fNrx =par[6]/kRad2Deg; // Nasmyth rotator displacement, horizontan + fNrx =par[6] /kRad2Deg; // Nasmyth rotator displacement, horizontan case 6: - fAw =par[5]/kRad2Deg; // Azimuth Axis Misalignment (E-W) + fAw =par[5] /kRad2Deg; // Azimuth Axis Misalignment (E-W) case 5: - fAn =par[4]/kRad2Deg; // Azimuth Axis Misalignment (N-S) + fAn =par[4] /kRad2Deg; // Azimuth Axis Misalignment (N-S) case 4: - fCa =par[3]/kRad2Deg; // Left-Right Collimation Error + fCa =par[3] /kRad2Deg; // Left-Right Collimation Error case 3: - fNpae =par[2]/kRad2Deg; // Az-El Nonperpendicularity + fNpae =par[2] /kRad2Deg; // Az-El Nonperpendicularity case 2: - fIe =par[1]/kRad2Deg; // Index Error in Elevation + fIe =par[1] /kRad2Deg; // Index Error in Elevation case 1: - fIa =par[0]/kRad2Deg; // Index Error in Azimuth + fIa =par[0] /kRad2Deg; // Index Error in Azimuth } } @@ -312,35 +314,35 @@ { case 16: - par[15]=fMagic2*kRad2Deg; // + par[15]=kRad2Deg*fMagic2; // case 15: - par[14]=fMagic1*kRad2Deg; // + par[14]=kRad2Deg*fMagic1; // case 14: - par[13]=fAcec*kRad2Deg; // Azimuth Centering Error (cos) + par[13]=kRad2Deg*fAcec; // Azimuth Centering Error (cos) case 13: - par[12]=fEcec*kRad2Deg; // Elevation Centering Error (cos) + par[12]=kRad2Deg*fEcec; // Elevation Centering Error (cos) case 12: - par[11]=fAces*kRad2Deg; // Azimuth Centering Error (sin) + par[11]=kRad2Deg*fAces; // Azimuth Centering Error (sin) case 11: - par[10]=fEces*kRad2Deg; // Elevation Centering Error (sin) + par[10]=kRad2Deg*fEces; // Elevation Centering Error (sin) case 10: - par[9]=fCry*kRad2Deg; // Alt/Az Coude Displacement (E-W) + par[9] =kRad2Deg*fCry; // Alt/Az Coude Displacement (E-W) case 9: - par[8]=fCrx*kRad2Deg; // Alt/Az Coude Displacement (N-S) + par[8] =kRad2Deg*fCrx; // Alt/Az Coude Displacement (N-S) case 8: - par[7]=fNry*kRad2Deg; // Nasmyth rotator displacement, vertical + par[7] =kRad2Deg*fNry; // Nasmyth rotator displacement, vertical case 7: - par[6]=fNrx*kRad2Deg; // Nasmyth rotator displacement, horizontan + par[6] =kRad2Deg*fNrx; // Nasmyth rotator displacement, horizontan case 6: - par[5]=fAw*kRad2Deg; // Azimuth Axis Misalignment (E-W) + par[5] =kRad2Deg*fAw; // Azimuth Axis Misalignment (E-W) case 5: - par[4]=fAn*kRad2Deg; // Azimuth Axis Misalignment (N-S) + par[4] =kRad2Deg*fAn; // Azimuth Axis Misalignment (N-S) case 4: - par[3]=fCa*kRad2Deg; // Left-Right Collimation Error + par[3] =kRad2Deg*fCa; // Left-Right Collimation Error case 3: - par[2]=fNpae*kRad2Deg; // Az-El Nonperpendicularity + par[2] =kRad2Deg*fNpae; // Az-El Nonperpendicularity case 2: - par[1]=fIe*kRad2Deg; // Index Error in Elevation + par[1] =kRad2Deg*fIe; // Index Error in Elevation case 1: - par[0]=fIa*kRad2Deg; // Index Error in Azimuth + par[0] =kRad2Deg*fIa; // Index Error in Azimuth } } @@ -349,5 +351,5 @@ { if (n<0) - n = m.GetNumPars(); + n = 16; Int_t ierflg = 0; @@ -357,50 +359,34 @@ case 16: m.mnparm(15,"MAGIC2", fMagic2*kRad2Deg, 1, -360, 360, ierflg); - // cout << "Init 3 CA: " << fCa << endl; case 15: m.mnparm(14,"MAGIC1", fMagic1*kRad2Deg, 1, -360, 360, ierflg); - // cout << "Init 3 CA: " << fCa << endl; case 14: m.mnparm(13,"ACEC", fAcec*kRad2Deg, 1, -360, 360, ierflg); - // cout << "Init 3 CA: " << fCa << endl; case 13: m.mnparm(12,"ECEC", fEcec*kRad2Deg, 1, -360, 360, ierflg); - // cout << "Init 3 CA: " << fCa << endl; case 12: m.mnparm(11,"ACES", fAcec*kRad2Deg, 1, -360, 360, ierflg); - // cout << "Init 3 CA: " << fCa << endl; case 11: m.mnparm(10,"ECES", fEcec*kRad2Deg, 1, -360, 360, ierflg); - // cout << "Init 3 CA: " << fCa << endl; case 10: m.mnparm(9, "CRY", fCry*kRad2Deg, 1, -360, 360, ierflg); - // cout << "Init 3 CA: " << fCa << endl; case 9: m.mnparm(8, "CRX", fCrx*kRad2Deg, 1, -360, 360, ierflg); - // cout << "Init 3 CA: " << fCa << endl; case 8: m.mnparm(7, "NRY", fNry*kRad2Deg, 1, -360, 360, ierflg); - // cout << "Init 3 CA: " << fCa << endl; case 7: m.mnparm(6, "NRX", fNrx*kRad2Deg, 1, -360, 360, ierflg); - // cout << "Init 3 CA: " << fCa << endl; case 6: m.mnparm(5, "AW", fAw*kRad2Deg, 1, -360, 360, ierflg); - // cout << "Init 3 CA: " << fCa << endl; case 5: m.mnparm(4, "AN", fAn*kRad2Deg, 1, -360, 360, ierflg); - // cout << "Init 3 CA: " << fCa << endl; case 4: m.mnparm(3, "CA", fCa*kRad2Deg, 1, -360, 360, ierflg); - // cout << "Init 3 CA: " << fCa << endl; case 3: m.mnparm(2, "NPAE", fNpae*kRad2Deg, 1, -360, 360, ierflg); - // cout << "Init 2 NPAE: " << fNpae << endl; case 2: - m.mnparm(1, "IE", fIe*kRad2Deg, 1, -360, 360, ierflg); - // cout << "Init 1 IE: " << fIe << endl; + m.mnparm(1, "IE", fIe*kRad2Deg, 1, -90, 90, ierflg); case 1: m.mnparm(0, "IA", fIa*kRad2Deg, 1, -360, 360, ierflg); - // cout << "Init 0 IA: " << fIa << endl; } } @@ -408,5 +394,5 @@ void MBending::GetMinuitParameters(TMinuit &m, Int_t n=-1) { - if (n<0) + if (n<0 || n>m.GetNumPars()) n = m.GetNumPars(); @@ -465,5 +451,24 @@ } } - +/* +void FormatPar(TMinuit &m, Int_t n) +{ + Double_t par, err; + m.GetParameter(n, par, err); + + int expp = (int)log10(par); + int expe = (int)log10(err); + + if (err<2*pow(10, expe)) + expe--; + + Int_t exp = expe>expp ? expp : expe; + + par = (int)(par/pow(10, exp)) * pow(10, exp); + err = (int)(err/pow(10, exp)) * pow(10, exp); + + cout << par << " +- " << err << flush; +} +*/ void MBending::PrintMinuitParameters(TMinuit &m, Int_t n=-1) const { @@ -471,4 +476,6 @@ n = m.GetNumPars(); + cout << setprecision(3); + Double_t par, err; @@ -477,50 +484,50 @@ case 16: m.GetParameter(15, par, err); - cout << " 15 MAGIC2: " << par << " +- " << err << endl; + cout << " 15 MAGIC2: " << setw(8) << par << " +- " << setw(4) << err << endl; case 15: m.GetParameter(14, par, err); - cout << " 14 MAGIC1: " << par << " +- " << err << endl; + cout << " 14 MAGIC1: " << setw(8) << par << " +- " << setw(4) << err << endl; case 14: m.GetParameter(13, par, err); - cout << " 13 ACEC: " << par << " +- " << err << endl; + cout << " 13 ACEC: " << setw(8) << par << " +- " << setw(4) << err << " Azimuth Centering Error (cos)" << endl; case 13: m.GetParameter(12, par, err); - cout << " 12 ECEC: " << par << " +- " << err << endl; + cout << " 12 ECEC: " << setw(8) << par << " +- " << setw(4) << err << " Elevation Centering Error (cos)" << endl; case 12: m.GetParameter(11, par, err); - cout << " 11 ACES: " << par << " +- " << err << endl; + cout << " 11 ACES: " << setw(8) << par << " +- " << setw(4) << err << " Azimuth Centering Error (sin)" << endl; case 11: m.GetParameter(10, par, err); - cout << " 10 ECES: " << par << " +- " << err << endl; + cout << " 10 ECES: " << setw(8) << par << " +- " << setw(4) << err << " Elevation Centering Error (sin)" << endl; case 10: m.GetParameter(9, par, err); - cout << " 9 CRY: " << par << " +- " << err << endl; + cout << " 9 CRY: " << setw(8) << par << " +- " << setw(4) << err << " Alt/Az Coude Displacement (E-W)" << endl; case 9: m.GetParameter(8, par, err); - cout << " 8 CRX: " << par << " +- " << err << endl; + cout << " 8 CRX: " << setw(8) << par << " +- " << setw(4) << err << " Alt/Az Coude Displacement (N-S)" << endl; case 8: m.GetParameter(7, par, err); - cout << " 7 NRY: " << par << " +- " << err << endl; + cout << " 7 NRY: " << setw(8) << par << " +- " << setw(4) << err << " Nasmyth rotator displacement, vertical" << endl; case 7: m.GetParameter(6, par, err); - cout << " 6 NRX: " << par << " +- " << err << endl; + cout << " 6 NRX: " << setw(8) << par << " +- " << setw(4) << err << " Nasmyth rotator displacement, horizontan" << endl; case 6: m.GetParameter(5, par, err); - cout << " 5 AW: " << par << " +- " << err << endl; + cout << " 5 AW: " << setw(8) << par << " +- " << setw(4) << err << " Azimuth Axis Misalignment (E-W)" << endl; case 5: m.GetParameter(4, par, err); - cout << " 4 AN: " << par << " +- " << err << endl; + cout << " 4 AN: " << setw(8) << par << " +- " << setw(4) << err << " Azimuth Axis Misalignment (N-S)" << endl; case 4: m.GetParameter(3, par, err); - cout << " 3 CA: " << par << " +- " << err << endl; + cout << " 3 CA: " << setw(8) << par << " +- " << setw(4) << err << " Left-Right Collimation Error" << endl; case 3: m.GetParameter(2, par, err); - cout << " 2 NPAE: " << par << " +- " << err << endl; + cout << " 2 NPAE: " << setw(8) << par << " +- " << setw(4) << err << " Az-El Nonperpendicularity" << endl; case 2: m.GetParameter(1, par, err); - cout << " 1 IE: " << par << " +- " << err << endl; + cout << " 1 IE: " << setw(8) << par << " +- " << setw(4) << err << " Index Error Elevation (Offset)" << endl; case 1: m.GetParameter(0, par, err); - cout << " 0 IA: " << par << " +- " << err << endl; - } -} + cout << " 0 IA: " << setw(8) << par << " +- " << setw(4) << err << " Index Error Azimuth (Offset)" << endl; + } +} Index: /trunk/MagicSoft/Cosy/tpoint/tpointfit.C =================================================================== --- /trunk/MagicSoft/Cosy/tpoint/tpointfit.C (revision 1809) +++ /trunk/MagicSoft/Cosy/tpoint/tpointfit.C (revision 1810) @@ -1,4 +1,5 @@ #include #include +#include #include @@ -21,12 +22,4 @@ // Sekans = 1/cos // - -void MyAdjust(AltAz &p, Double_t *par) -{ - // p[rad] - MBending bend; - bend.SetParameters(par); // [deg] - p=bend(p); -} class Set : public TObject @@ -54,16 +47,13 @@ Double_t GetResidual() const { - Double_t daz = fRawEl-fStarEl; - Double_t del = fRawAz-fStarAz; - - //Double_t rzd = TMath::Pi()/2-fRawEl; - Double_t rzd = TMath::Pi()/2-fOrigStarEl; + 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(2*rzd+del)*sin2; - - Double_t dist = acos(d); + Double_t d = cos(del)*cos2 - cos(2*fOrigStarEl+del)*sin2; + + Double_t dist = acos(d); return dist * 180 / TMath::Pi(); @@ -79,4 +69,5 @@ 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(); } @@ -102,12 +93,9 @@ fStarAz = p.Az(); } - - void Adjust(const MBending &bend, void (*fcn)(AltAz &zdaz, Double_t *par)) - { - AltAz star = GetStarAltAz(); - - AltAz p = bend(star, fcn); - fStarEl = p.Alt(); - fStarAz = p.Az(); + void AdjustBack(const MBending &bend) + { + AltAz p = bend.CorrectBack(GetRawAltAz()); + fRawEl = p.Alt(); + fRawAz = p.Az(); } }; @@ -116,8 +104,8 @@ { Double_t v[4]; + fin >> v[0]; + fin >> v[1]; fin >> v[2]; fin >> v[3]; - fin >> v[0]; - fin >> v[1]; set.fStarAz = v[0]*TMath::Pi()/180; @@ -140,5 +128,5 @@ MBending bend; - bend.SetParameters(par/*, npar*/); // Set Parameters [deg] to MBending + bend.SetParameters(par); // Set Parameters [deg] to MBending for (int i=0; i1000) + scale = -1; + if (scale>0) { Double_t d = r1-r0; - if (d<0) d=-d; - r0 -= scale*d; - r1 += scale*d; + r0 += scale*d; + r1 -= scale*d; d = phi1-phi0; - if (d<0) d=-d; - phi0 -= scale*d; - phi1 += scale*d; + phi0 += scale*d; + phi1 -= scale*d; } @@ -290,19 +296,7 @@ gCoordinates.SetOwner(); - TH2F hdaz("dAz", "\\Delta Az vs. Alt", 32, 0, 90, 32, -1, 1); - TH2F hdzd("dZd", "\\Delta Alt vs. Az", 32, 0, 360, 32, -1, 1); - TH1F hres1("Res1", " Residuals before correction ", 10, 0, 180); TH1F hres2("Res2", " Residuals after correction ", 10, 0, 1); - TH2F h2res1("Res2D1", " Arb. Residuals before correction ", 32, 0, 360, 10, 0, 90); - TH2F h2res2("Res2D2", " Arb. Residuals after correction ", 32, 0, 360, 10, 0, 90); - - hdaz.SetXTitle("Zd [\\circ]"); - hdaz.SetYTitle("\\Delta Az [\\circ]"); - - hdzd.SetXTitle("Az [\\circ]"); - hdzd.SetYTitle("\\Delta Zd [\\circ]"); - hres1.SetXTitle("\\Delta [\\circ]"); hres1.SetYTitle("Counts"); @@ -313,6 +307,14 @@ TGraph gdaz; TGraph gdzd; - - ifstream fin("tpoint/tpoint2.txt"); + TGraph gaz; + TGraph gzd; + + gdaz.SetTitle(" \\Delta Az vs. Zd "); + gdzd.SetTitle(" \\Delta Zd vs. Az "); + + gaz.SetTitle(" \\Delta Az vs. Az "); + gzd.SetTitle(" \\Delta Zd vs. Zd "); + + ifstream fin("tpoint/tpoint3.txt"); while (fin && fin.get()!='\n'); @@ -346,65 +348,116 @@ bending.SetMinuitParameters(minuit, 16); // Init Parameters [deg] - Int_t ierflg = 0; - - // minuit.FixParameter(2); // NPAE - // minuit.FixParameter(3); // CA - // minuit.FixParameter(4); // AN - // minuit.FixParameter(5); // AW - // minuit.FixParameter(6); // NRX - // minuit.FixParameter(7); // NRY + //minuit.FixParameter(0); //(1) IA + //minuit.FixParameter(1); //(1) IE + minuit.FixParameter(2); //(2) NPAE + //minuit.FixParameter(3); //(1) CA + minuit.FixParameter(4); //(2) AN + //minuit.FixParameter(5); //(1) AW + + minuit.FixParameter(6); // NRX + minuit.FixParameter(7); // NRY minuit.FixParameter(8); // CRX minuit.FixParameter(9); // CRY - // minuit.FixParameter(10); // ECES - // minuit.FixParameter(11); // ACES - // minuit.FixParameter(12); // ECEC - // minuit.FixParameter(13); // ACEC - // minuit.FixParameter(14); // MAGIC1 + minuit.FixParameter(10); // ECES + minuit.FixParameter(11); //(2) ACES + minuit.FixParameter(12); // ECEC + minuit.FixParameter(13); //(2) ACEC + minuit.FixParameter(14); // MAGIC1 minuit.FixParameter(15); // MAGIC2 + //minuit.Command("SHOW PARAMETERS"); + //minuit.Command("SHOW LIMITS"); + cout << endl; + + Int_t ierflg = 0; ierflg = minuit.Migrad(); - cout << endl << "Migrad returns " << ierflg << endl << endl; - minuit.Command("SHOW LIMITS"); - - cout << endl; - - /* - minuit.FixParameter(0); - minuit.FixParameter(1); - minuit.Release(2); - - ierflg = minuit.Migrad(); - cout << endl << "Migrad returns " << ierflg << endl << endl; - */ - + cout << "Migrad returns " << ierflg << endl; + // minuit.Release(2); + ierflg = minuit.Migrad(); + cout << "Migrad returns " << ierflg << endl << endl; + + // + // Get Fit Results + // bending.GetMinuitParameters(minuit); bending.PrintMinuitParameters(minuit); bending.Save("bending_magic.txt"); + + // + // Make a copy of all list entries + // TList list; list.SetOwner(); - for (int i=0; i180) + dz -= 360; static int j=0; - gdzd.SetPoint(j, set0.GetStarAz(), set0.GetDEl()); - gdaz.SetPoint(j, 90-set0.GetStarEl(), set0.GetDAz()); + gdzd.SetPoint(j, za.Az()/* set0.GetStarAz()*/, set0.GetDZd()); + gaz.SetPoint( j, za.Az()/* set0.GetStarAz()*/, dz); + gdaz.SetPoint(j, za.Zd()/*90-set0.GetStarEl()*/, dz); + gzd.SetPoint( j, za.Zd()/*90-set0.GetStarEl()*/, set0.GetDZd()); j++; } + // + // Check for overflows + // + const Stat_t ov = hres2.GetBinContent(hres2.GetNbinsX()+1); + if (ov>0) + cout << "WARNING: " << ov << " overflows in residuals." << endl; + + // + // 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; iDivide(3,2); + c1->Divide(2,2); + + c1->cd(2); + gdaz.DrawClone("A*"); c1->cd(1); - hdaz.DrawCopy("cont"); - gdaz.DrawClone("*"); + gaz.DrawClone("A*"); + + c1->cd(3); + gdzd.DrawClone("A*"); c1->cd(4); - hdzd.DrawCopy("cont"); - gdzd.DrawClone("*"); + gzd.DrawClone("A*"); + + + + c1=new TCanvas("c2", "c", 800, 800); + c1->Divide(2,2); c1->cd(2); @@ -429,25 +493,27 @@ hres1.DrawCopy(); - c1->cd(5); + c1->cd(4); hres2.SetLineColor(kBlue); hres2.DrawCopy(); -// c1->cd(3); -// gPad->SetTheta(90); -// gPad->SetPhi(90); -// h2res1.DrawCopy("surf1pol"); -// gPad->Modified(); -// gPad->Update(); -// for (int i=0; icd(6); -// gPad->SetTheta(90); -// gPad->SetPhi(90); -// h2res1.DrawCopy("surf1pol"); -// gPad->Modified(); -// gPad->Update(); -// for (int i=0; icd(1); + gPad->SetTheta(90); + gPad->SetPhi(90); + TH2F h2res1("Res2D1", " Arb. Residuals before correction (scaled) ", 32, 0, 360, 10, 0, 90); + h2res1.DrawCopy("surf1pol"); + gPad->Modified(); + gPad->Update(); + for (int i=0; icd(3); + 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