Changeset 1810
- Timestamp:
- 03/11/03 13:52:48 (22 years ago)
- Location:
- trunk/MagicSoft/Cosy
- Files:
-
- 13 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Cosy/Changelog
r1803 r1810 1 1 -*-*- END -*-*- 2 2003/03/11 - Daniela Dorner, Thomas Bretz: 3 4 * base/MStar.h: 5 - added Compare 6 7 * base/MStarList.cc: 8 - some small bugfixes 9 10 * base/MStarList.h: 11 - added Sort 12 - added Expand 13 14 * base/timer.[h,cc]: 15 - Added GetTimeval 16 17 * gui/MGAccuracy.cc: 18 - Exchanged zd/az in calculation of Residual!!! 19 20 * gui/MGPngReader.[h,cc]: 21 - set default lim mag to 7.0 22 - added new ouput for the pointing position fPZdAz 23 - added/fixed TrackingError/CalcTrackingError 24 - changed Filter2 to CaosFilter 25 - reordered starguider stuff in Execute 26 - changed color of circles 27 28 * main/MBending.[h,cc]: 29 - removed MAGIC1 and MAGIC2 30 - removed '-' from writing 31 - fixed some bugs in the enumerations of the coefficients 32 - added some formating option for output 33 34 * tpoint/tpointfit.C: 35 - removed usage of MyAdjust 36 - fixed the Calculation of the residuals 37 - fixed reading 38 - added some correction in case of an overflow (360deg/0deg) 39 - fixed drawing 40 - added second Migrad turn... 41 - changed the screen and graphical output 42 43 * videodev/CaosFilter.[h,cc]: 44 - changed RemoveTwins to accept a radius 45 46 2 47 3 48 2003/03/02 - Daniela Dorner, Thomas Bretz (LaPalma): -
trunk/MagicSoft/Cosy/base/MStar.h
r1691 r1810 21 21 void Set(Double_t mx, Double_t my) { fX=mx; fY=my; } 22 22 23 Int_t Compare(const TObject *obj) const 24 { 25 const MStar *const s = (MStar*)obj; 26 27 if (fMag<s->fMag) 28 return -1; 29 30 if (fMag>s->fMag) 31 return 1; 32 33 return 0; 34 } 35 36 Bool_t IsSortable() const { return kTRUE; } 37 38 23 39 ClassDef(MStar, 1) 24 40 }; -
trunk/MagicSoft/Cosy/base/MStarList.cc
r1691 r1810 1 1 #include "MStarList.h" 2 3 #include <iostream.h> 2 4 3 5 void MStarList::RemoveTwins(Double_t radius) … … 18 20 return; 19 21 22 fStars.RemoveAt(idx); 23 20 24 MStarListIter Next(this, *first, radius); 21 Delete(idx);22 25 23 26 MStar *pos; … … 29 32 mx += pos->GetX(); 30 33 my += pos->GetY(); 31 Delete(pos);34 fStars.Remove(pos); 32 35 cnt++; 33 36 } -
trunk/MagicSoft/Cosy/base/MStarList.h
r1760 r1810 38 38 39 39 void RemoveTwins(Double_t radius); 40 41 void Sort() { fStars.Sort(); } 42 43 void Expand(int n) { fStars.Expand(n); } 40 44 }; 41 45 -
trunk/MagicSoft/Cosy/base/timer.cc
r1793 r1810 26 26 // fDiv = fmod((fMs+fSecs)/(60*60*24), 1.0); 27 27 } 28 29 void Timer::GetTimeval(struct timeval *tv) const 30 { 31 tv->tv_sec = fSecs; 32 tv->tv_usec = fMs; 33 } 34 28 35 /* 29 36 void Timer::Set(const long mjd) -
trunk/MagicSoft/Cosy/base/timer.h
r1760 r1810 33 33 void SetTimer(const struct timeval *tv); 34 34 35 void GetTimeval(struct timeval *tv) const; 36 35 37 int GetSecs() { return fSecs; } 36 38 -
trunk/MagicSoft/Cosy/gui/MGAccuracy.cc
r1760 r1810 195 195 aaz *= d2r; 196 196 197 const double el = TMath::Pi()/2-pzd; 198 197 199 const double dphi2 = aaz/2.; 198 200 const double cos2 = cos(dphi2)*cos(dphi2); 199 201 const double sin2 = sin(dphi2)*sin(dphi2); 200 const double d = cos(azd)*cos2 - cos(2* pzd+azd)*sin2;202 const double d = cos(azd)*cos2 - cos(2*el+azd)*sin2; 201 203 202 204 double dist = acos(d); … … 250 252 251 253 UpdateCross(x, y); 252 UpdateText(pos.Zd(), acc. Az(), acc.Zd());254 UpdateText(pos.Zd(), acc.Zd(), acc.Az()); 253 255 254 256 SetModified(); -
trunk/MagicSoft/Cosy/gui/MGPngReader.cc
r1802 r1810 3 3 #include <fstream.h> // ifstream 4 4 #include <iostream.h> // cout 5 #include <iomanip.h> // cout 5 6 6 7 #include <TTimer.h> 7 8 8 9 #include <TGMenu.h> 10 #include <TGLabel.h> 9 11 #include <TSystem.h> 10 12 #include <TGSplitter.h> // TGHorizontal3DLine … … 13 15 #include "MGImage.h" 14 16 #include "MGCoordinates.h" 17 #include "MGDispCoordinates.h" 15 18 16 19 #include "coord.h" … … 255 258 fList->Add(fLimMag); 256 259 257 fSao->SetLimitMag( 8.0);260 fSao->SetLimitMag(7.0); 258 261 259 262 fInterpol = new MGPopupMenu(p); … … 314 317 fList->Add(fCRaDec); 315 318 316 fCZdAz = new MGCoordinates(this, kETypeZdAz );319 fCZdAz = new MGCoordinates(this, kETypeZdAz, kFALSE); 317 320 fCZdAz->Move(240+12+10, fMenu->GetDefaultHeight()+584); 318 321 AddFrame(fCZdAz); 319 322 fList->Add(fCZdAz); 320 323 324 fPZdAz = new MGCoordinates(this, kETypeZdAz, kFALSE); 325 fPZdAz->Move(240+12+10, fMenu->GetDefaultHeight()+630); 326 AddFrame(fPZdAz); 327 fList->Add(fPZdAz); 328 329 TGLabel *l = new TGLabel(this, "Arb.-Sky Pos"); 330 l->SetTextJustify(kTextLeft); 331 l->Move(480+32, fMenu->GetDefaultHeight()+590); 332 AddFrame(l); 333 fList->Add(l); 334 335 l = new TGLabel(this, "arcsec/pix"); 336 l->SetTextJustify(kTextLeft); 337 l->Move(605, fMenu->GetDefaultHeight()+619); 338 AddFrame(l); 339 fList->Add(l); 340 341 l = new TGLabel(this, "Pointing Pos"); 342 l->SetTextJustify(kTextLeft); 343 l->Move(480+32, fMenu->GetDefaultHeight()+655); 344 AddFrame(l); 345 fList->Add(l); 346 321 347 const Double_t pixsize = 23.4; 322 348 … … 328 354 fPixSize = new TGTextEntry(this, txt, IDM_kPixSize); 329 355 fPixSize->SetAlignment(kTextCenterX); 330 fPixSize->Move( 600, fMenu->GetDefaultHeight()+584);356 fPixSize->Move(547, fMenu->GetDefaultHeight()+617); 331 357 AddFrame(fPixSize); 332 358 fList->Add(fPixSize); … … 357 383 MapSubwindows(); 358 384 MapWindow(); 385 386 //------------------------------------------------------------ 387 // XY xy(3.819444, 24.05333); 388 // fCRaDec->SetCoordinates(xy); 389 // fRaDec->Set(xy.X()*360/24, xy.Y()); 390 //------------------------------------------------------------ 359 391 } 360 392 361 393 MGPngReader::MGPngReader(MObservatory::LocationName_t obs) 362 : TGMainFrame(gClient->GetRoot(), 768, 7 00), fFile(NULL), fDx((768-kZOOM)/2), fDy((512-kZOOM)/2)394 : TGMainFrame(gClient->GetRoot(), 768, 740), fFile(NULL), fDx((768-kZOOM)/2), fDy((512-kZOOM)/2) 363 395 { 364 396 fSao = new StarCatalog(obs); … … 811 843 } 812 844 813 void MGPngReader::CalcTrackingError(MStarList &spots, MStarList &stars) 814 { 845 ZdAz MGPngReader::TrackingError(TArrayF &x, TArrayF &y, TArrayF &mag) const 846 { 847 // 848 // Viewable area (FIXME: AZ) 849 // 850 TH2F h("Hist", "dX/dY", 77, -768/2-.5, 768/2+.5, 58, -576/2-.5, 576/2+.5); // 3 851 852 /* 853 TH1F hmag("HistMag", "Mag", 19, 0, 100); 854 for (int i=0; i<mag.GetSize(); i++) 855 hmag.Fill(mag[i]); 856 */ 857 858 // 859 // Search for matching Magnitudes 860 // 861 for (int i=0; i<mag.GetSize(); i++) 862 { 863 if (mag[i]>48-15 && mag[i]<48+15) 864 h.Fill(x[i], y[i]); 865 } 866 867 // 868 // Serach for an excess in the histogram 869 // 870 Int_t mx, my, dummy; 871 h.GetMaximumBin(mx, my, dummy); 872 873 const double xmax = h.GetXaxis()->GetBinCenter(mx); 874 const double dx = h.GetXaxis()->GetBinWidth(mx); 875 876 const double ymax = h.GetYaxis()->GetBinCenter(my); 877 const double dy = h.GetYaxis()->GetBinWidth(my); 878 879 cout << setprecision(3); 880 cout << "Cut-XY: " << xmax << " +- " << dx << " / " << ymax << " +- " << dy << endl; 881 882 TGraph g; 883 for (int i=0; i<mag.GetSize(); i++) 884 { 885 if (!(x[i]>xmax-dx && x[i]<xmax+dx && 886 y[i]>ymax-dy && y[i]<ymax+dy && 887 mag[i]>48-15 && mag[i]<48+15)) 888 continue; 889 890 g.SetPoint(g.GetN(), x[i], y[i]); 891 } 892 893 cout << "Offset-XY: " << g.GetMean(1) << " +- " << g.GetRMS(1) << " / "; 894 cout << g.GetMean(2) << " +- " << g.GetRMS(2) << endl; 895 896 AltAz pos0 = fSao->CalcAltAzFromPix(768/2, 576/2)*kRad2Deg; 897 AltAz pos1 = fSao->CalcAltAzFromPix(768/2+g.GetMean(1), 576/2+g.GetMean(2))*kRad2Deg; 898 899 pos1 -= pos0; 900 901 ofstream fout("tracking_error.txt"); 902 fout << setprecision(10) << fSao->GetMjd()-52000 << " " << -pos1.Alt() << " " << pos1.Az() << endl; 903 fout.close(); 904 905 return ZdAz(-pos1.Alt(), pos1.Az()); 906 } 907 908 void MGPngReader::CalcTrackingError(Leds &leds, MStarList &stars) 909 { 910 const Int_t max = leds.GetEntries(); 911 815 912 if (stars.GetRealEntries() < 3) 816 913 { … … 819 916 } 820 917 821 if ( spots.GetRealEntries()< 1)918 if (max < 1) 822 919 { 823 920 cout << "Sorry, less than 1 detected spot in FOV!" << endl; … … 825 922 } 826 923 827 Int_t idx = 0; 828 829 MStarList sortedspots; 924 stars.Sort(); // Sort by magnitude 925 926 TString str = "data/tracking_"; 927 str += fSao->GetMjd()-52000; 928 str += ".txt"; 929 930 ofstream fout(str); 931 932 TArrayF x, y, mag; 933 934 Int_t num = 0; 935 936 // FIXME: Is predifined value 3 a good idea? 830 937 831 938 MStar *star; 832 MStar *spot;833 939 MStarListIter NextStar(&stars); 834 MStarListIter NextSpot(&spots); 835 836 while ((spot=NextSpot())) 837 { 838 AltAz aa = fSao->CalcAltAzFromPix(spot->GetX(), spot->GetY()); 839 spot->Set(aa.Az(), aa.Alt()); 840 } 841 842 while ((star=NextStar())) 843 { 844 AltAz aa = fSao->CalcAltAzFromPix(star->GetX(), star->GetY()); 845 star->Set(aa.Az(), aa.Alt()); 846 847 const double aaz = star->GetX(); 848 const double dphi2 = aaz/2.; 849 const double cos2 = cos(dphi2)*cos(dphi2); 850 const double sin2 = sin(dphi2)*sin(dphi2); 851 852 Double_t min = 800; 853 854 NextSpot.Reset(); 855 while ((spot=NextSpot())) 940 while ((star=NextStar()) && num++<max+3) 941 { 942 TIter NextSp(&leds); 943 Led *spot=NULL; 944 while ((spot=(Led*)NextSp())) 856 945 { 857 const double pzd = TMath::Pi()/2-spot->GetY(); 858 const double azd = TMath::Pi()/2-star->GetY(); 859 860 const double d = cos(azd)*cos2 - cos(2*pzd+azd)*sin2; 861 862 const Double_t dist = acos(d); 863 864 if (dist>=min) 865 continue; 866 867 min = dist; 868 sortedspots.AddAt(idx, spot->GetX(), spot->GetY(), spot->GetMag()); 946 const XY dpos(spot->GetX()-star->GetX(), spot->GetY()-star->GetY()); 947 948 const Int_t idx = x.GetSize(); 949 950 x.Set(idx+1); 951 y.Set(idx+1); 952 mag.Set(idx+1); 953 954 x.AddAt(dpos.X(), idx); 955 y.AddAt(dpos.Y(), idx); 956 mag.AddAt(spot->GetMag()/star->GetMag(), idx); 957 958 if (fout) 959 fout << x[idx] << " " << y[idx] << " " << mag[idx] << endl; 869 960 } 870 if (min>768) 871 { 872 cout << "ERROR!!!!!!!!" << endl; 873 return; 874 } 875 idx++; 876 } 877 878 // 879 // Now we have in sortedspots the entries with the shortest distances 880 // to the corresponding ones in stars. 881 // Now calculate the tracking error. 882 // 883 NextStar.Reset(); 884 MStarListIter NextSpot2(&sortedspots); 885 886 Double_t meanx=0; 887 Double_t meany=0; 888 889 while ((star=NextStar())) 890 { 891 spot = NextSpot2(); 892 893 meanx += star->GetX() - spot->GetX(); 894 meany += star->GetY() - spot->GetY(); 895 } 896 897 meanx /= idx; 898 meany /= idx; 899 900 cout << "Tracking Error: dAlt=" << meany*180/TMath::Pi(); 901 cout << "° dAz=" << meanx*180/TMath::Pi() << "° (calculated"; 902 cout << " with " << idx << " stars/spots)" << endl; 903 } 904 961 } 962 963 ZdAz d = TrackingError(x, y, mag); 964 965 // 966 // Calculated offsets 967 // 968 969 // round= floor(x+.5) 970 cout << "Offset-ZdAz: " << d.Zd()*60 << "' / " << d.Az()*60 << "'" << endl; 971 cout << "Offset-ZdAz: " << d.Zd()/360*16384 << " / " << d.Az()/360*16384 << " (SE) " << endl; 972 973 // 974 // Current Pointing position 975 // 976 ZdAz cpos = fSao->GetZdAz()-d; 977 fPZdAz->SetCoordinates(cpos); 978 } 905 979 906 980 … … 945 1019 946 1020 MStarList spots; 1021 1022 /* 947 1023 if (fDisplay->IsEntryChecked(IDM_kStarguider)) 948 1024 Filter2::Execute(spots, c); 949 1025 else 950 if (fDisplay->IsEntryChecked(IDM_kFilter)) 951 Filter::Execute(c); 1026 */ 1027 if (fDisplay->IsEntryChecked(IDM_kFilter)) 1028 Filter::Execute(c); 952 1029 953 1030 if (fDisplay->IsEntryChecked(IDM_kCaosFilter)) … … 1031 1108 if (fDisplay->IsEntryChecked(IDM_kCatalog)) 1032 1109 { 1110 Timer time(tm); 1111 1112 GetCoordinates(); 1113 1114 MStarList stars; 1115 fSao->GetStars(stars, time.GetMjd(), *fRaDec); 1116 1117 if (fDisplay->IsEntryChecked(IDM_kStarguider)) 1118 { 1119 Leds leds; 1120 CaosFilter::Execute(c, leds, 1); 1121 1122 cout << "Found: " << leds.GetEntries() << " leds" << endl; 1123 1124 CaosFilter::RemoveTwins(leds, 3); 1125 leds.Compress(); 1126 1127 cout << "Rest: " << leds.GetEntries() << " leds" << endl; 1128 1129 CalcTrackingError(leds, stars); 1130 } 1131 1132 byte cimg[768*576]; 1133 fSao->GetImg(c, cimg, stars); 1134 1033 1135 DrawCircle(c, 0.5); 1034 1136 DrawCircle(c, 1.0); 1035 1137 DrawCircle(c, 1.5); 1036 1138 1037 byte cimg[768*576]; 1038 1039 GetCoordinates(); 1040 1041 Timer time(tm); 1042 1043 MStarList stars; 1044 fSao->GetStars(stars, time.GetMjd(), *fRaDec); 1045 fSao->GetImg(c, cimg, stars); 1046 //fSao->GetImg(c, cimg, time.CalcMjd(), *fRaDec); 1139 fCZdAz->SetCoordinates(fSao->GetZdAz()); 1047 1140 1048 1141 fImage->DrawColImg(c, cimg); 1049 1050 fCZdAz->SetCoordinates(fSao->GetZdAz());1051 1052 if (fDisplay->IsEntryChecked(IDM_kStarguider))1053 CalcTrackingError(spots, stars);1054 1142 } 1055 1143 else … … 1067 1155 { 1068 1156 const int dy = (int)sqrt(rpix*rpix-dx*dx); 1069 img[cx+dx + (cy-dy)*768] = 0x d0;1070 img[cx+dx + (cy+dy)*768] = 0x d0;1071 img[cx-dy + (cy+dx)*768] = 0x d0;1072 img[cx+dy + (cy+dx)*768] = 0x d0;1157 img[cx+dx + (cy-dy)*768] = 0x40; 1158 img[cx+dx + (cy+dy)*768] = 0x40; 1159 img[cx-dy + (cy+dx)*768] = 0x40; 1160 img[cx+dy + (cy+dx)*768] = 0x40; 1073 1161 } 1074 1162 } -
trunk/MagicSoft/Cosy/gui/MGPngReader.h
r1802 r1810 13 13 #include "MGImage.h" 14 14 15 #include <TH1.h> 16 #include <TH2.h> 17 #include <TGraph.h> 18 #include <TCanvas.h> 15 #include "coord.h" 19 16 20 class AltAz; 21 class RaDec; 17 class TArrayF; 18 class TH1F; 19 class TH2F; 20 class TGraph; 22 21 23 22 class TTimer; … … 87 86 MGCoordinates *fCZdAz; 88 87 88 MGCoordinates *fPZdAz; 89 89 90 TGTextEntry *fPixSize; 90 91 … … 104 105 void Toggle(MGPopupMenu *p, UInt_t id); 105 106 void GetCoordinates(); 106 void CalcTrackingError(MStarList &, MStarList &); 107 void CalcTrackingError(Leds &, MStarList &); 108 ZdAz TrackingError(TArrayF &alt, TArrayF &az, TArrayF &mag) const; 107 109 108 110 void InitHists(); -
trunk/MagicSoft/Cosy/main/MBending.cc
r1808 r1810 107 107 if (str=="ECEC") fEcec = val; 108 108 if (str=="ACEC") fAcec = val; 109 if (str=="MAGIC1") fMagic1 = val;110 if (str=="MAGIC2") fMagic2 = val;109 //if (str=="MAGIC1") fMagic1 = val; 110 //if (str=="MAGIC2") fMagic2 = val; 111 111 112 112 fin >> val; … … 143 143 fout << "S 00 000000 000000 0000000" << endl; 144 144 fout << setprecision(8); 145 fout << " IA " << -kRad2Deg*fIa << " -1" << endl;146 fout << " IE " << -kRad2Deg*fIe << " -1" << endl;147 fout << " CA " << -kRad2Deg*fNpae << " -1" << endl;148 fout << " NPAE " << -kRad2Deg*fCa << " -1" << endl;149 fout << " AN " << -kRad2Deg*fAn << " -1" << endl;150 fout << " AW " << -kRad2Deg*fAw << " -1" << endl;151 fout << " NRX " << -kRad2Deg*fNrx << " -1" << endl;152 fout << " NRY " << -kRad2Deg*fNry << " -1" << endl;153 fout << " CRX " << -kRad2Deg*fCrx << " -1" << endl;154 fout << " CRY " << -kRad2Deg*fCry << " -1" << endl;155 fout << " ECES " << -kRad2Deg*fEces << " -1" << endl;156 fout << " ACES " << -kRad2Deg*fAces << " -1" << endl;157 fout << " ECEC " << -kRad2Deg*fEcec << " -1" << endl;158 fout << " ACEC " << -kRad2Deg*fAcec << " -1" << endl;159 fout << " MAGIC1 " << -kRad2Deg*fMagic1 << " -1" << endl;160 fout << " MAGIC2 " << -kRad2Deg*fMagic2 << " -1" << endl;145 fout << " IA " << kRad2Deg*fIa << " -1" << endl; 146 fout << " IE " << kRad2Deg*fIe << " -1" << endl; 147 fout << " CA " << kRad2Deg*fNpae << " -1" << endl; 148 fout << " NPAE " << kRad2Deg*fCa << " -1" << endl; 149 fout << " AN " << kRad2Deg*fAn << " -1" << endl; 150 fout << " AW " << kRad2Deg*fAw << " -1" << endl; 151 fout << " NRX " << kRad2Deg*fNrx << " -1" << endl; 152 fout << " NRY " << kRad2Deg*fNry << " -1" << endl; 153 fout << " CRX " << kRad2Deg*fCrx << " -1" << endl; 154 fout << " CRY " << kRad2Deg*fCry << " -1" << endl; 155 fout << " ECES " << kRad2Deg*fEces << " -1" << endl; 156 fout << " ACES " << kRad2Deg*fAces << " -1" << endl; 157 fout << " ECEC " << kRad2Deg*fEcec << " -1" << endl; 158 fout << " ACEC " << kRad2Deg*fAcec << " -1" << endl; 159 //fout << " MAGIC1 " << -kRad2Deg*fMagic1 << " -1" << endl; 160 //fout << " MAGIC2 " << -kRad2Deg*fMagic2 << " -1" << endl; 161 161 fout << "END" << endl; 162 162 } … … 173 173 p += CEC; 174 174 175 // const AltAz MAGIC(0, -fMagic1*tan(p.Alt())-fMagic2); 176 // p += MAGIC; 177 175 178 const AltAz CRX(-fCrx*sin(p.Az()-p.Alt()), fCrx*cos(p.Az()-p.Alt())/cos(p.Alt())); 176 179 const AltAz CRY(-fCry*cos(p.Az()-p.Alt()), -fCry*sin(p.Az()-p.Alt())/cos(p.Alt())); … … 183 186 p += NRY; 184 187 185 const AltAz MAGIC(-fMagic1*sin(p.Az()-fMagic2), 0);186 p += MAGIC;188 // const AltAz MAGIC(-fMagic1*sin(p.Az()-fMagic2), 0); 189 // p += MAGIC; 187 190 188 191 const AltAz AW( fAw*sin(p.Az()), -fAw*cos(p.Az())*tan(p.Alt())); … … 191 194 p += AN; 192 195 193 // const AltAz MAGIC(-fMagic1*sin(p.Az()-fMagic2), 0); 194 // p += MAGIC; 195 196 // const AltAz MAGIC(-fMagic1*sin(p.Az()-fMagic2), 0); 197 // p += MAGIC; 196 198 197 199 const AltAz CA(0, -fCa/cos(p.Alt())); … … 227 229 p -= AW; 228 230 229 const AltAz MAGIC(-fMagic1*sin(p.Az()-fMagic2), 0);230 p -= MAGIC;231 // const AltAz MAGIC(-fMagic1*sin(p.Az()-fMagic2), 0); 232 // p -= MAGIC; 231 233 232 234 const AltAz NRY(fNry*cos(p.Alt()), -fNry*tan(p.Alt())); … … 285 287 fEces =par[10]/kRad2Deg; // Elevation Centering Error (sin) 286 288 case 10: 287 fCry =par[9] /kRad2Deg;// Alt/Az Coude Displacement (E-W)289 fCry =par[9] /kRad2Deg; // Alt/Az Coude Displacement (E-W) 288 290 case 9: 289 fCrx =par[8] /kRad2Deg;// Alt/Az Coude Displacement (N-S)291 fCrx =par[8] /kRad2Deg; // Alt/Az Coude Displacement (N-S) 290 292 case 8: 291 fNry =par[7] /kRad2Deg;// Nasmyth rotator displacement, vertical293 fNry =par[7] /kRad2Deg; // Nasmyth rotator displacement, vertical 292 294 case 7: 293 fNrx =par[6] /kRad2Deg;// Nasmyth rotator displacement, horizontan295 fNrx =par[6] /kRad2Deg; // Nasmyth rotator displacement, horizontan 294 296 case 6: 295 fAw =par[5] /kRad2Deg;// Azimuth Axis Misalignment (E-W)297 fAw =par[5] /kRad2Deg; // Azimuth Axis Misalignment (E-W) 296 298 case 5: 297 fAn =par[4] /kRad2Deg;// Azimuth Axis Misalignment (N-S)299 fAn =par[4] /kRad2Deg; // Azimuth Axis Misalignment (N-S) 298 300 case 4: 299 fCa =par[3] /kRad2Deg;// Left-Right Collimation Error301 fCa =par[3] /kRad2Deg; // Left-Right Collimation Error 300 302 case 3: 301 fNpae =par[2] /kRad2Deg;// Az-El Nonperpendicularity303 fNpae =par[2] /kRad2Deg; // Az-El Nonperpendicularity 302 304 case 2: 303 fIe =par[1] /kRad2Deg;// Index Error in Elevation305 fIe =par[1] /kRad2Deg; // Index Error in Elevation 304 306 case 1: 305 fIa =par[0] /kRad2Deg;// Index Error in Azimuth307 fIa =par[0] /kRad2Deg; // Index Error in Azimuth 306 308 } 307 309 } … … 312 314 { 313 315 case 16: 314 par[15]= fMagic2*kRad2Deg; //316 par[15]=kRad2Deg*fMagic2; // 315 317 case 15: 316 par[14]= fMagic1*kRad2Deg; //318 par[14]=kRad2Deg*fMagic1; // 317 319 case 14: 318 par[13]= fAcec*kRad2Deg; // Azimuth Centering Error (cos)320 par[13]=kRad2Deg*fAcec; // Azimuth Centering Error (cos) 319 321 case 13: 320 par[12]= fEcec*kRad2Deg; // Elevation Centering Error (cos)322 par[12]=kRad2Deg*fEcec; // Elevation Centering Error (cos) 321 323 case 12: 322 par[11]= fAces*kRad2Deg; // Azimuth Centering Error (sin)324 par[11]=kRad2Deg*fAces; // Azimuth Centering Error (sin) 323 325 case 11: 324 par[10]= fEces*kRad2Deg; // Elevation Centering Error (sin)326 par[10]=kRad2Deg*fEces; // Elevation Centering Error (sin) 325 327 case 10: 326 par[9] =fCry*kRad2Deg;// Alt/Az Coude Displacement (E-W)328 par[9] =kRad2Deg*fCry; // Alt/Az Coude Displacement (E-W) 327 329 case 9: 328 par[8] =fCrx*kRad2Deg;// Alt/Az Coude Displacement (N-S)330 par[8] =kRad2Deg*fCrx; // Alt/Az Coude Displacement (N-S) 329 331 case 8: 330 par[7] =fNry*kRad2Deg;// Nasmyth rotator displacement, vertical332 par[7] =kRad2Deg*fNry; // Nasmyth rotator displacement, vertical 331 333 case 7: 332 par[6] =fNrx*kRad2Deg;// Nasmyth rotator displacement, horizontan334 par[6] =kRad2Deg*fNrx; // Nasmyth rotator displacement, horizontan 333 335 case 6: 334 par[5] =fAw*kRad2Deg;// Azimuth Axis Misalignment (E-W)336 par[5] =kRad2Deg*fAw; // Azimuth Axis Misalignment (E-W) 335 337 case 5: 336 par[4] =fAn*kRad2Deg;// Azimuth Axis Misalignment (N-S)338 par[4] =kRad2Deg*fAn; // Azimuth Axis Misalignment (N-S) 337 339 case 4: 338 par[3] =fCa*kRad2Deg;// Left-Right Collimation Error340 par[3] =kRad2Deg*fCa; // Left-Right Collimation Error 339 341 case 3: 340 par[2] =fNpae*kRad2Deg; // Az-El Nonperpendicularity342 par[2] =kRad2Deg*fNpae; // Az-El Nonperpendicularity 341 343 case 2: 342 par[1] =fIe*kRad2Deg;// Index Error in Elevation344 par[1] =kRad2Deg*fIe; // Index Error in Elevation 343 345 case 1: 344 par[0] =fIa*kRad2Deg;// Index Error in Azimuth346 par[0] =kRad2Deg*fIa; // Index Error in Azimuth 345 347 } 346 348 } … … 349 351 { 350 352 if (n<0) 351 n = m.GetNumPars();353 n = 16; 352 354 353 355 Int_t ierflg = 0; … … 357 359 case 16: 358 360 m.mnparm(15,"MAGIC2", fMagic2*kRad2Deg, 1, -360, 360, ierflg); 359 // cout << "Init 3 CA: " << fCa << endl;360 361 case 15: 361 362 m.mnparm(14,"MAGIC1", fMagic1*kRad2Deg, 1, -360, 360, ierflg); 362 // cout << "Init 3 CA: " << fCa << endl;363 363 case 14: 364 364 m.mnparm(13,"ACEC", fAcec*kRad2Deg, 1, -360, 360, ierflg); 365 // cout << "Init 3 CA: " << fCa << endl;366 365 case 13: 367 366 m.mnparm(12,"ECEC", fEcec*kRad2Deg, 1, -360, 360, ierflg); 368 // cout << "Init 3 CA: " << fCa << endl;369 367 case 12: 370 368 m.mnparm(11,"ACES", fAcec*kRad2Deg, 1, -360, 360, ierflg); 371 // cout << "Init 3 CA: " << fCa << endl;372 369 case 11: 373 370 m.mnparm(10,"ECES", fEcec*kRad2Deg, 1, -360, 360, ierflg); 374 // cout << "Init 3 CA: " << fCa << endl;375 371 case 10: 376 372 m.mnparm(9, "CRY", fCry*kRad2Deg, 1, -360, 360, ierflg); 377 // cout << "Init 3 CA: " << fCa << endl;378 373 case 9: 379 374 m.mnparm(8, "CRX", fCrx*kRad2Deg, 1, -360, 360, ierflg); 380 // cout << "Init 3 CA: " << fCa << endl;381 375 case 8: 382 376 m.mnparm(7, "NRY", fNry*kRad2Deg, 1, -360, 360, ierflg); 383 // cout << "Init 3 CA: " << fCa << endl;384 377 case 7: 385 378 m.mnparm(6, "NRX", fNrx*kRad2Deg, 1, -360, 360, ierflg); 386 // cout << "Init 3 CA: " << fCa << endl;387 379 case 6: 388 380 m.mnparm(5, "AW", fAw*kRad2Deg, 1, -360, 360, ierflg); 389 // cout << "Init 3 CA: " << fCa << endl;390 381 case 5: 391 382 m.mnparm(4, "AN", fAn*kRad2Deg, 1, -360, 360, ierflg); 392 // cout << "Init 3 CA: " << fCa << endl;393 383 case 4: 394 384 m.mnparm(3, "CA", fCa*kRad2Deg, 1, -360, 360, ierflg); 395 // cout << "Init 3 CA: " << fCa << endl;396 385 case 3: 397 386 m.mnparm(2, "NPAE", fNpae*kRad2Deg, 1, -360, 360, ierflg); 398 // cout << "Init 2 NPAE: " << fNpae << endl;399 387 case 2: 400 m.mnparm(1, "IE", fIe*kRad2Deg, 1, -360, 360, ierflg); 401 // cout << "Init 1 IE: " << fIe << endl; 388 m.mnparm(1, "IE", fIe*kRad2Deg, 1, -90, 90, ierflg); 402 389 case 1: 403 390 m.mnparm(0, "IA", fIa*kRad2Deg, 1, -360, 360, ierflg); 404 // cout << "Init 0 IA: " << fIa << endl;405 391 } 406 392 } … … 408 394 void MBending::GetMinuitParameters(TMinuit &m, Int_t n=-1) 409 395 { 410 if (n<0 )396 if (n<0 || n>m.GetNumPars()) 411 397 n = m.GetNumPars(); 412 398 … … 465 451 } 466 452 } 467 453 /* 454 void FormatPar(TMinuit &m, Int_t n) 455 { 456 Double_t par, err; 457 m.GetParameter(n, par, err); 458 459 int expp = (int)log10(par); 460 int expe = (int)log10(err); 461 462 if (err<2*pow(10, expe)) 463 expe--; 464 465 Int_t exp = expe>expp ? expp : expe; 466 467 par = (int)(par/pow(10, exp)) * pow(10, exp); 468 err = (int)(err/pow(10, exp)) * pow(10, exp); 469 470 cout << par << " +- " << err << flush; 471 } 472 */ 468 473 void MBending::PrintMinuitParameters(TMinuit &m, Int_t n=-1) const 469 474 { … … 471 476 n = m.GetNumPars(); 472 477 478 cout << setprecision(3); 479 473 480 Double_t par, err; 474 481 … … 477 484 case 16: 478 485 m.GetParameter(15, par, err); 479 cout << " 15 MAGIC2: " << par << " +- " <<err << endl;486 cout << " 15 MAGIC2: " << setw(8) << par << " +- " << setw(4) << err << endl; 480 487 case 15: 481 488 m.GetParameter(14, par, err); 482 cout << " 14 MAGIC1: " << par << " +- " <<err << endl;489 cout << " 14 MAGIC1: " << setw(8) << par << " +- " << setw(4) << err << endl; 483 490 case 14: 484 491 m.GetParameter(13, par, err); 485 cout << " 13 ACEC: " << par << " +- " << err<< endl;492 cout << " 13 ACEC: " << setw(8) << par << " +- " << setw(4) << err << " Azimuth Centering Error (cos)" << endl; 486 493 case 13: 487 494 m.GetParameter(12, par, err); 488 cout << " 12 ECEC: " << par << " +- " << err<< endl;495 cout << " 12 ECEC: " << setw(8) << par << " +- " << setw(4) << err << " Elevation Centering Error (cos)" << endl; 489 496 case 12: 490 497 m.GetParameter(11, par, err); 491 cout << " 11 ACES: " << par << " +- " << err<< endl;498 cout << " 11 ACES: " << setw(8) << par << " +- " << setw(4) << err << " Azimuth Centering Error (sin)" << endl; 492 499 case 11: 493 500 m.GetParameter(10, par, err); 494 cout << " 10 ECES: " << par << " +- " << err<< endl;501 cout << " 10 ECES: " << setw(8) << par << " +- " << setw(4) << err << " Elevation Centering Error (sin)" << endl; 495 502 case 10: 496 503 m.GetParameter(9, par, err); 497 cout << " 9 CRY: " << par << " +- " << err<< endl;504 cout << " 9 CRY: " << setw(8) << par << " +- " << setw(4) << err << " Alt/Az Coude Displacement (E-W)" << endl; 498 505 case 9: 499 506 m.GetParameter(8, par, err); 500 cout << " 8 CRX: " << par << " +- " << err<< endl;507 cout << " 8 CRX: " << setw(8) << par << " +- " << setw(4) << err << " Alt/Az Coude Displacement (N-S)" << endl; 501 508 case 8: 502 509 m.GetParameter(7, par, err); 503 cout << " 7 NRY: " << par << " +- " << err<< endl;510 cout << " 7 NRY: " << setw(8) << par << " +- " << setw(4) << err << " Nasmyth rotator displacement, vertical" << endl; 504 511 case 7: 505 512 m.GetParameter(6, par, err); 506 cout << " 6 NRX: " << par << " +- " << err<< endl;513 cout << " 6 NRX: " << setw(8) << par << " +- " << setw(4) << err << " Nasmyth rotator displacement, horizontan" << endl; 507 514 case 6: 508 515 m.GetParameter(5, par, err); 509 cout << " 5 AW: " << par << " +- " << err<< endl;516 cout << " 5 AW: " << setw(8) << par << " +- " << setw(4) << err << " Azimuth Axis Misalignment (E-W)" << endl; 510 517 case 5: 511 518 m.GetParameter(4, par, err); 512 cout << " 4 AN: " << par << " +- " << err<< endl;519 cout << " 4 AN: " << setw(8) << par << " +- " << setw(4) << err << " Azimuth Axis Misalignment (N-S)" << endl; 513 520 case 4: 514 521 m.GetParameter(3, par, err); 515 cout << " 3 CA: " << par << " +- " << err<< endl;522 cout << " 3 CA: " << setw(8) << par << " +- " << setw(4) << err << " Left-Right Collimation Error" << endl; 516 523 case 3: 517 524 m.GetParameter(2, par, err); 518 cout << " 2 NPAE: " << par << " +- " << err<< endl;525 cout << " 2 NPAE: " << setw(8) << par << " +- " << setw(4) << err << " Az-El Nonperpendicularity" << endl; 519 526 case 2: 520 527 m.GetParameter(1, par, err); 521 cout << " 1 IE: " << par << " +- " << err<< endl;528 cout << " 1 IE: " << setw(8) << par << " +- " << setw(4) << err << " Index Error Elevation (Offset)" << endl; 522 529 case 1: 523 530 m.GetParameter(0, par, err); 524 cout << " 0 IA: " << par << " +- " << err<< endl;525 } 526 } 531 cout << " 0 IA: " << setw(8) << par << " +- " << setw(4) << err << " Index Error Azimuth (Offset)" << endl; 532 } 533 } -
trunk/MagicSoft/Cosy/tpoint/tpointfit.C
r1806 r1810 1 1 #include <fstream.h> 2 2 #include <iostream.h> 3 #include <iomanip.h> 3 4 4 5 #include <TF1.h> … … 21 22 // Sekans = 1/cos 22 23 // 23 24 void MyAdjust(AltAz &p, Double_t *par)25 {26 // p[rad]27 MBending bend;28 bend.SetParameters(par); // [deg]29 p=bend(p);30 }31 24 32 25 class Set : public TObject … … 54 47 Double_t GetResidual() const 55 48 { 56 Double_t daz = fRawEl-fStarEl; 57 Double_t del = fRawAz-fStarAz; 58 59 //Double_t rzd = TMath::Pi()/2-fRawEl; 60 Double_t rzd = TMath::Pi()/2-fOrigStarEl; 49 Double_t del = fRawEl-fStarEl; 50 Double_t daz = fRawAz-fStarAz; 61 51 62 52 Double_t dphi2 = daz/2.; 63 53 Double_t cos2 = cos(dphi2)*cos(dphi2); 64 54 Double_t sin2 = sin(dphi2)*sin(dphi2); 65 Double_t d = cos(del)*cos2 - cos(2*rzd+del)*sin2;66 67 Double_t dist = acos(d);55 Double_t d = cos(del)*cos2 - cos(2*fOrigStarEl+del)*sin2; 56 57 Double_t dist = acos(d); 68 58 69 59 return dist * 180 / TMath::Pi(); … … 79 69 80 70 Double_t GetDEl() const { return (fRawEl-fStarEl)*180/TMath::Pi(); } 71 Double_t GetDZd() const { return -GetDEl(); } 81 72 Double_t GetDAz() const { return (fRawAz-fStarAz)*180/TMath::Pi(); } 82 73 Double_t GetStarEl() const { return fStarEl*180/TMath::Pi(); } … … 102 93 fStarAz = p.Az(); 103 94 } 104 105 void Adjust(const MBending &bend, void (*fcn)(AltAz &zdaz, Double_t *par)) 106 { 107 AltAz star = GetStarAltAz(); 108 109 AltAz p = bend(star, fcn); 110 fStarEl = p.Alt(); 111 fStarAz = p.Az(); 95 void AdjustBack(const MBending &bend) 96 { 97 AltAz p = bend.CorrectBack(GetRawAltAz()); 98 fRawEl = p.Alt(); 99 fRawAz = p.Az(); 112 100 } 113 101 }; … … 116 104 { 117 105 Double_t v[4]; 106 fin >> v[0]; 107 fin >> v[1]; 118 108 fin >> v[2]; 119 109 fin >> v[3]; 120 fin >> v[0];121 fin >> v[1];122 110 123 111 set.fStarAz = v[0]*TMath::Pi()/180; … … 140 128 141 129 MBending bend; 142 bend.SetParameters(par /*, npar*/); // Set Parameters [deg] to MBending130 bend.SetParameters(par); // Set Parameters [deg] to MBending 143 131 144 132 for (int i=0; i<gCoordinates.GetSize(); i++) … … 146 134 Set set = *(Set*)gCoordinates.At(i); 147 135 148 //set.Adjust(bend); 149 set.Adjust(bend, MyAdjust); 136 set.Adjust(bend); 150 137 151 138 Double_t err = 1; … … 200 187 return; 201 188 } 202 189 /* 203 190 if (r0<0) 204 191 { … … 212 199 } 213 200 214 /* 215 phi0 = fmod(phi0+360, 360); 216 phi1 = fmod(phi1+360, 360); 201 phi0 = fmod(phi0+360, 360); 202 phi1 = fmod(phi1+360, 360); 203 204 if (phi1-phi0<-180) 205 phi1+=360; 217 206 */ 218 219 207 TLine line; 220 208 line.SetLineWidth(2); … … 268 256 Double_t phi1 = set.GetStarAz(); 269 257 258 if (r0<0) 259 { 260 r0 = -r0; 261 phi0 += 180; 262 } 263 if (r1<0) 264 { 265 r1 = -r1; 266 phi1 += 180; 267 } 268 269 phi0 = fmod(phi0+360, 360); 270 phi1 = fmod(phi1+360, 360); 271 272 if (phi1-phi0<-180) 273 phi1+=360; 274 275 if (scale<0 || scale>1000) 276 scale = -1; 277 270 278 if (scale>0) 271 279 { 272 280 Double_t d = r1-r0; 273 if (d<0) d=-d; 274 r0 -= scale*d; 275 r1 += scale*d; 281 r0 += scale*d; 282 r1 -= scale*d; 276 283 d = phi1-phi0; 277 if (d<0) d=-d; 278 phi0 -= scale*d; 279 phi1 += scale*d; 284 phi0 += scale*d; 285 phi1 -= scale*d; 280 286 } 281 287 … … 290 296 gCoordinates.SetOwner(); 291 297 292 TH2F hdaz("dAz", "\\Delta Az vs. Alt", 32, 0, 90, 32, -1, 1);293 TH2F hdzd("dZd", "\\Delta Alt vs. Az", 32, 0, 360, 32, -1, 1);294 295 298 TH1F hres1("Res1", " Residuals before correction ", 10, 0, 180); 296 299 TH1F hres2("Res2", " Residuals after correction ", 10, 0, 1); 297 300 298 TH2F h2res1("Res2D1", " Arb. Residuals before correction ", 32, 0, 360, 10, 0, 90);299 TH2F h2res2("Res2D2", " Arb. Residuals after correction ", 32, 0, 360, 10, 0, 90);300 301 hdaz.SetXTitle("Zd [\\circ]");302 hdaz.SetYTitle("\\Delta Az [\\circ]");303 304 hdzd.SetXTitle("Az [\\circ]");305 hdzd.SetYTitle("\\Delta Zd [\\circ]");306 307 301 hres1.SetXTitle("\\Delta [\\circ]"); 308 302 hres1.SetYTitle("Counts"); … … 313 307 TGraph gdaz; 314 308 TGraph gdzd; 315 316 ifstream fin("tpoint/tpoint2.txt"); 309 TGraph gaz; 310 TGraph gzd; 311 312 gdaz.SetTitle(" \\Delta Az vs. Zd "); 313 gdzd.SetTitle(" \\Delta Zd vs. Az "); 314 315 gaz.SetTitle(" \\Delta Az vs. Az "); 316 gzd.SetTitle(" \\Delta Zd vs. Zd "); 317 318 ifstream fin("tpoint/tpoint3.txt"); 317 319 318 320 while (fin && fin.get()!='\n'); … … 346 348 bending.SetMinuitParameters(minuit, 16); // Init Parameters [deg] 347 349 348 Int_t ierflg = 0; 349 350 // minuit.FixParameter(2); // NPAE 351 // minuit.FixParameter(3); // CA 352 // minuit.FixParameter(4); // AN 353 // minuit.FixParameter(5); // AW 354 // minuit.FixParameter(6); // NRX 355 // minuit.FixParameter(7); // NRY 350 //minuit.FixParameter(0); //(1) IA 351 //minuit.FixParameter(1); //(1) IE 352 minuit.FixParameter(2); //(2) NPAE 353 //minuit.FixParameter(3); //(1) CA 354 minuit.FixParameter(4); //(2) AN 355 //minuit.FixParameter(5); //(1) AW 356 357 minuit.FixParameter(6); // NRX 358 minuit.FixParameter(7); // NRY 356 359 minuit.FixParameter(8); // CRX 357 360 minuit.FixParameter(9); // CRY 358 //minuit.FixParameter(10); // ECES359 // minuit.FixParameter(11); //ACES360 //minuit.FixParameter(12); // ECEC361 // minuit.FixParameter(13); //ACEC362 //minuit.FixParameter(14); // MAGIC1361 minuit.FixParameter(10); // ECES 362 minuit.FixParameter(11); //(2) ACES 363 minuit.FixParameter(12); // ECEC 364 minuit.FixParameter(13); //(2) ACEC 365 minuit.FixParameter(14); // MAGIC1 363 366 minuit.FixParameter(15); // MAGIC2 364 367 368 365 369 //minuit.Command("SHOW PARAMETERS"); 370 //minuit.Command("SHOW LIMITS"); 371 cout << endl; 372 373 Int_t ierflg = 0; 366 374 ierflg = minuit.Migrad(); 367 cout << endl << "Migrad returns " << ierflg << endl << endl; 368 minuit.Command("SHOW LIMITS"); 369 370 cout << endl; 371 372 /* 373 minuit.FixParameter(0); 374 minuit.FixParameter(1); 375 minuit.Release(2); 376 377 ierflg = minuit.Migrad(); 378 cout << endl << "Migrad returns " << ierflg << endl << endl; 379 */ 380 375 cout << "Migrad returns " << ierflg << endl; 376 // minuit.Release(2); 377 ierflg = minuit.Migrad(); 378 cout << "Migrad returns " << ierflg << endl << endl; 379 380 // 381 // Get Fit Results 382 // 381 383 bending.GetMinuitParameters(minuit); 382 384 bending.PrintMinuitParameters(minuit); 383 385 bending.Save("bending_magic.txt"); 384 386 387 388 // 389 // Make a copy of all list entries 390 // 385 391 TList list; 386 392 list.SetOwner(); 387 388 393 for (int i=0; i<gCoordinates.GetSize(); i++) 394 list.Add(new Set(*(Set*)gCoordinates.At(i))); 395 396 // 397 // Calculate correction and residuals 398 // 399 for (int i=0; i<gCoordinates.GetSize(); i++) 389 400 { 390 401 Set &set0 = *(Set*)gCoordinates.At(i); 391 402 392 list.Add(new Set(set0));403 ZdAz za = set0.GetStarZdAz()*kRad2Deg; 393 404 394 405 hres1.Fill(set0.GetResidual()); 395 //set0.Adjust(bending); 396 set0.Adjust(bending, MyAdjust); 406 set0.Adjust(bending); 397 407 hres2.Fill(set0.GetResidual()); 398 408 399 hdzd.Fill( set0.GetStarAz(), set0.GetDEl()); 400 hdaz.Fill(90-set0.GetStarEl(), set0.GetDAz()); 409 Double_t dz = fmod(set0.GetDAz()+720, 360); 410 if (dz>180) 411 dz -= 360; 401 412 402 413 static int j=0; 403 gdzd.SetPoint(j, set0.GetStarAz(), set0.GetDEl()); 404 gdaz.SetPoint(j, 90-set0.GetStarEl(), set0.GetDAz()); 414 gdzd.SetPoint(j, za.Az()/* set0.GetStarAz()*/, set0.GetDZd()); 415 gaz.SetPoint( j, za.Az()/* set0.GetStarAz()*/, dz); 416 gdaz.SetPoint(j, za.Zd()/*90-set0.GetStarEl()*/, dz); 417 gzd.SetPoint( j, za.Zd()/*90-set0.GetStarEl()*/, set0.GetDZd()); 405 418 j++; 406 419 } 407 420 421 // 422 // Check for overflows 423 // 424 const Stat_t ov = hres2.GetBinContent(hres2.GetNbinsX()+1); 425 if (ov>0) 426 cout << "WARNING: " << ov << " overflows in residuals." << endl; 427 428 // 429 // Print all data sets for which the backward correction is 430 // twice times worse than the residual gotten from the 431 // bending correction itself 432 // 408 433 cout << endl; 434 cout << "Checking backward correction (raw-->star):" << endl; 435 for (int i=0; i<gCoordinates.GetSize(); i++) 436 { 437 Set set0(*(Set*)list.At(i)); 438 Set set1(*(Set*)list.At(i)); 439 440 set0.AdjustBack(bending); 441 set1.Adjust(bending); 442 443 Double_t res0 = set0.GetResidual(); 444 Double_t res1 = set1.GetResidual(); 445 446 Double_t diff = fabs(res0-res1); 447 448 if (diff<hres2.GetMean()*0.66) 449 continue; 450 451 cout << "DBack: " << setw(7) << set0.GetStarZd() << " " << setw(8) << set0.GetStarAz() << ": "; 452 cout << "ResB="<< setw(7) << res0*60 << " ResF=" << setw(7) << res1*60 << " |ResB-ResF|=" << setw(7) << diff*60 << " arcmin" << endl; 453 } 454 cout << "OK." << endl; 455 456 // 457 // Print out the residual before and after correction in several 458 // units 459 // 460 cout << endl; 461 cout << gCoordinates.GetSize() << " data sets." << endl << endl; 409 462 cout << "Spead before: " << hres1.GetMean() << "deg +- " << hres1.GetRMS() << "deg" << endl; 410 463 cout << "Spead after: " << hres2.GetMean() << "deg +- " << hres2.GetRMS() << "deg" << endl; … … 414 467 cout << endl; 415 468 469 470 416 471 TCanvas *c1=new TCanvas("c1", "c"); 417 c1->Divide(3,2); 472 c1->Divide(2,2); 473 474 c1->cd(2); 475 gdaz.DrawClone("A*"); 418 476 419 477 c1->cd(1); 420 hdaz.DrawCopy("cont"); 421 gdaz.DrawClone("*"); 478 gaz.DrawClone("A*"); 479 480 c1->cd(3); 481 gdzd.DrawClone("A*"); 422 482 423 483 c1->cd(4); 424 hdzd.DrawCopy("cont"); 425 gdzd.DrawClone("*"); 484 gzd.DrawClone("A*"); 485 486 487 488 c1=new TCanvas("c2", "c", 800, 800); 489 c1->Divide(2,2); 426 490 427 491 c1->cd(2); … … 429 493 hres1.DrawCopy(); 430 494 431 c1->cd( 5);495 c1->cd(4); 432 496 hres2.SetLineColor(kBlue); 433 497 hres2.DrawCopy(); 434 498 435 // c1->cd(3); 436 // gPad->SetTheta(90); 437 // gPad->SetPhi(90); 438 // h2res1.DrawCopy("surf1pol"); 439 // gPad->Modified(); 440 // gPad->Update(); 441 // for (int i=0; i<gCoordinates.GetSize(); i++) 442 // DrawSet(gPad, *(Set*)list.At(i), 10./hres1.GetMean()); 443 444 // c1->cd(6); 445 // gPad->SetTheta(90); 446 // gPad->SetPhi(90); 447 // h2res1.DrawCopy("surf1pol"); 448 // gPad->Modified(); 449 // gPad->Update(); 450 // for (int i=0; i<gCoordinates.GetSize(); i++) 451 // DrawSet(gPad, *(Set*)gCoordinates.At(i), 10./hres2.GetMean()); 499 c1->cd(1); 500 gPad->SetTheta(90); 501 gPad->SetPhi(90); 502 TH2F h2res1("Res2D1", " Arb. Residuals before correction (scaled) ", 32, 0, 360, 10, 0, 90); 503 h2res1.DrawCopy("surf1pol"); 504 gPad->Modified(); 505 gPad->Update(); 506 for (int i=0; i<gCoordinates.GetSize(); i++) 507 DrawSet(gPad, *(Set*)list.At(i), 10./hres1.GetMean()); 508 509 c1->cd(3); 510 gPad->SetTheta(90); 511 gPad->SetPhi(90); 512 h2res1.SetTitle(" Arb. Residuals after correction (scaled) "); 513 h2res1.DrawCopy("surf1pol"); 514 gPad->Modified(); 515 gPad->Update(); 516 for (int i=0; i<gCoordinates.GetSize(); i++) 517 DrawSet(gPad, *(Set*)gCoordinates.At(i), 10./hres2.GetMean()); 452 518 453 519 } -
trunk/MagicSoft/Cosy/videodev/CaosFilter.cc
r1802 r1810 345 345 } 346 346 347 RemoveTwins(leds, 5); 348 } 349 350 void CaosFilter::RemoveTwins(Leds &leds, Double_t radius) 351 { 347 352 TIter Next1(&leds); 348 353 … … 368 373 const Double_t dy = y2-y1; 369 374 370 if (dx*dx+dy*dy< 5*5)375 if (dx*dx+dy*dy<radius*radius) 371 376 { 372 377 // FIXME: Interpolation -
trunk/MagicSoft/Cosy/videodev/CaosFilter.h
r1803 r1810 49 49 static void Execute(byte *img, Leds &leds, Double_t conv); 50 50 51 static void FilterLeds(Leds &leds); 51 static void FilterLeds(Leds &leds); 52 static void RemoveTwins(Leds &leds, Double_t radius); 52 53 53 54 ClassDef(CaosFilter, 0)
Note:
See TracChangeset
for help on using the changeset viewer.