Index: trunk/MagicSoft/Cosy/Makefile =================================================================== --- trunk/MagicSoft/Cosy/Makefile (revision 4106) +++ trunk/MagicSoft/Cosy/Makefile (revision 4107) @@ -145,4 +145,6 @@ @ln -sf ../../Mars/mbase/MGMenu.h base/MGMenu.h @ln -sf ../../Mars/mbase/MGMenu.cc base/MGMenu.cc + @ln -sf ../../Mars/mbase/MString.h base/MString.h + @ln -sf ../../Mars/mbase/MString.cc base/MString.cc @ln -sf ../../Mars/mgeom/MGeomMirror.h base/MGeomMirror.h @ln -sf ../../Mars/mgeom/MGeomMirror.cc base/MGeomMirror.cc @@ -217,8 +219,8 @@ @cp stars.txt /home/stesy/Cosy/ @chmod a-w /home/stesy/Cosy/stars.txt - @echo Copy SAO catalog - @chmod u+w /home/stesy/Cosy/sao/sao-* - @cp sao/sao-* /home/stesy/Cosy/sao/ - @chmod a-w /home/stesy/Cosy/sao/sao-* + @echo Copy PPM catalog + @chmod u+w /home/stesy/Cosy/ppm9.bin + @cp ppm9.bin /home/stesy/Cosy/ + @chmod a-w /home/stesy/Cosy/ppm9.bin @echo finished... Index: trunk/MagicSoft/Cosy/caos/Led.cc =================================================================== --- trunk/MagicSoft/Cosy/caos/Led.cc (revision 4106) +++ trunk/MagicSoft/Cosy/caos/Led.cc (revision 4107) @@ -7,4 +7,5 @@ #include "Ring.h" +#include "MString.h" ClassImp(Led); @@ -17,8 +18,9 @@ void Led::Print(Option_t *o=NULL) const { + MString str; cout << "Led: "; - cout << "x=" << Form("%5.1f", fX) << "+-" << fDx << ", "; - cout << "y=" << Form("%5.1f", fY) << "+-" << fDy << ", "; - cout << "phi=" << Form("%6.1f", fPhi) << "+-" << fDphi << ", "; + cout << "x=" << str.Print("%5.1f", fX) << "+-" << fDx << ", "; + cout << "y=" << str.Print("%5.1f", fY) << "+-" << fDy << ", "; + cout << "phi=" << str.Print("%6.1f", fPhi) << "+-" << fDphi << ", "; cout << "mag=" << fMag << endl; } Index: trunk/MagicSoft/Cosy/caos/Makefile =================================================================== --- trunk/MagicSoft/Cosy/caos/Makefile (revision 4106) +++ trunk/MagicSoft/Cosy/caos/Makefile (revision 4107) @@ -20,5 +20,6 @@ # @endcode -INCLUDES = -I. -I.. +INCLUDES = -I. -I.. -I../base +#base: MString # @code Index: trunk/MagicSoft/Cosy/caos/Ring.cc =================================================================== --- trunk/MagicSoft/Cosy/caos/Ring.cc (revision 4106) +++ trunk/MagicSoft/Cosy/caos/Ring.cc (revision 4107) @@ -7,9 +7,10 @@ #include "Rings.h" +#include "MString.h" ClassImp(Ring); -Ring::Ring() : - fX(0), fY(0), fR(0), fPhi(0), fDx(-1), fDy(-1), fDr(-1), fDphi(-1) +Ring::Ring(Double_t x, Double_t y) : + fX(x), fY(y), fR(0), fPhi(0), fDx(-1), fDy(-1), fDr(-1), fDphi(-1) { } @@ -114,8 +115,9 @@ void Ring::Print(Option_t *o=NULL) const { + MString str; cout << "Ring: "; - cout << "x=" << Form("%5.1f", fX) << "+-" << Form("%.1f", fDx) << ", "; - cout << "y=" << Form("%5.1f", fY) << "+-" << Form("%.1f", fDy) << ", "; - cout << "r=" << Form("%5.1f", fR) << "+-" << Form("%.1f", fDr) << ", "; + cout << "x=" << str.Print("%5.1f", fX) << "+-" << Form("%.1f", fDx) << ", "; + cout << "y=" << str.Print("%5.1f", fY) << "+-" << Form("%.1f", fDy) << ", "; + cout << "r=" << str.Print("%5.1f", fR) << "+-" << Form("%.1f", fDr) << ", "; cout << "phi=" << fPhi << "+-" << fDphi << endl; } Index: trunk/MagicSoft/Cosy/caos/Ring.h =================================================================== --- trunk/MagicSoft/Cosy/caos/Ring.h (revision 4106) +++ trunk/MagicSoft/Cosy/caos/Ring.h (revision 4107) @@ -32,5 +32,5 @@ public: - Ring(); + Ring(Double_t x=0, Double_t y=0); Double_t GetX() const { return fX; } Index: trunk/MagicSoft/Cosy/main/MCosy.cc =================================================================== --- trunk/MagicSoft/Cosy/main/MCosy.cc (revision 4106) +++ trunk/MagicSoft/Cosy/main/MCosy.cc (revision 4107) @@ -883,5 +883,5 @@ cout << " * Y: " << kGear.Y() << "*" << kResRE.Y()/4 << "/" << kResSE.Y() << "=4*" << kGearTot.Y() << "/" << kResSE.Y() << endl; } - +/* void MCosy::InitSync() { @@ -899,5 +899,5 @@ fMac3->StartPosSync(); } - +*/ void MCosy::TalkThreadSeTest() { @@ -1042,5 +1042,5 @@ } - InitSync(); + //InitSync(); /*** FOR DEMO MODE ***/ Index: trunk/MagicSoft/Cosy/main/MTracking.cc =================================================================== --- trunk/MagicSoft/Cosy/main/MTracking.cc (revision 4106) +++ trunk/MagicSoft/Cosy/main/MTracking.cc (revision 4107) @@ -550,7 +550,10 @@ // time.Zd(fCosy->fMac2->GetMjd()); - Double_t mjd1 = fCosy->fZd1->GetMjd(); - Double_t mjd2 = fCosy->fZd2->GetMjd(); - Double_t mjd0 = fCosy->fAz->GetMjd(); + //Double_t mjd1 = fCosy->fZd1->GetMjd(); + //Double_t mjd2 = fCosy->fZd2->GetMjd(); + //Double_t mjd0 = fCosy->fAz->GetMjd(); + + Double_t mjdaz = fCosy->fMac1->GetPdoMjd();//mjd0; + Double_t mjdzd = fCosy->fMac2->GetPdoMjd();//TMath::Max(mjd1, mjd2); // get current position of shaftencoders (interpolated @@ -565,5 +568,5 @@ if (phca1 || phca2) { - ZdAz dummy = sla.CalcZdAz(fCosy->fRaDec, fCosy->fMac2->GetPdoMjd());//TMath::Max(mjd1, mjd2)); + ZdAz dummy = sla.CalcZdAz(fCosy->fRaDec, mjdzd); dummy = fCosy->AlignTrackingPos(dummy); dummy = fCosy->fBending(dummy); @@ -572,6 +575,5 @@ if (phcaz) { - const Double_t mjd = fCosy->fMac1->GetMjd(); - ZdAz dummy = sla.CalcZdAz(fCosy->fRaDec, fCosy->fMac1->GetPdoMjd());//mjd0); + ZdAz dummy = sla.CalcZdAz(fCosy->fRaDec, mjdaz); dummy = fCosy->AlignTrackingPos(dummy); dummy = fCosy->fBending(dummy); Index: trunk/MagicSoft/Cosy/tpoint/gui.C =================================================================== --- trunk/MagicSoft/Cosy/tpoint/gui.C (revision 4106) +++ trunk/MagicSoft/Cosy/tpoint/gui.C (revision 4107) @@ -1,5 +1,9 @@ +#include +#include #include #include #include + +#include #include @@ -38,23 +42,55 @@ 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(); + fStarAz(saz*TMath::DegToRad()), + fStarEl(sel*TMath::DegToRad()), + fRawAz(raz*TMath::DegToRad()), + fRawEl(rel*TMath::DegToRad()) + { + } + + Double_t GetResidual(Double_t *err=0) const + { + /* + TVector3 v1, v2; + v1.SetMagThetaPhi(1, TMath::Pi()/2-fRawEl, fRawAz); + v2.SetMagThetaPhi(1, TMath::Pi()/2-fStarEl, fStarAz); + + return v1.Angle(v2)*TMath::RadToDeg(); + */ + + const Double_t del = fRawEl-fStarEl; + const Double_t daz = fRawAz-fStarAz; + /* + const Double_t dphi2 = daz/2.; + const Double_t cos2 = cos(dphi2)*cos(dphi2); + const Double_t sin2 = sin(dphi2)*sin(dphi2); + const Double_t d = cos(del)*cos2 - cos(fRawEl+fStarEl)*sin2; + */ + + const Double_t d = cos(del) - cos(fRawEl)*cos(fStarEl)*(1.-cos(daz)); + + if (err) + { + // Error of one pixel in the CCD + const Double_t e1 = 32./3600*TMath::DegToRad() * 0.5; + + // Error of one SE unit + const Double_t e2 = 360./16384*TMath::DegToRad() * 0.5; + + const Double_t e11 = sin(del)+cos(fRawEl)*sin(fStarEl)*(1-cos(daz)); + const Double_t e12 = cos(fRawEl)*cos(fStarEl)*sin(daz); + + const Double_t e21 = -sin(del)+sin(fRawEl)*cos(fStarEl)*(1-cos(daz)); + const Double_t e22 = -cos(fRawEl)*cos(fStarEl)*sin(daz); + + const Double_t err1 = sqrt(1-d*d); + const Double_t err2 = (e11*e11 + e12*e12)*e1*e1; + const Double_t err3 = (e21*e21 + e22*e22)*e2*e2; + + *err = sqrt(err2+err3)/err1 * TMath::RadToDeg(); + } + + const Double_t dist = acos(d); + return dist * TMath::RadToDeg(); } @@ -67,13 +103,13 @@ } - Double_t GetDEl() const { return (fRawEl-fStarEl)*180/TMath::Pi(); } + Double_t GetDEl() const { return (fRawEl-fStarEl)*TMath::RadToDeg(); } 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(); } + Double_t GetDAz() const { return (fRawAz-fStarAz)*TMath::RadToDeg(); } + Double_t GetStarEl() const { return fStarEl*TMath::RadToDeg(); } + Double_t GetStarZd() const { return 90.-fStarEl*TMath::RadToDeg(); } + Double_t GetStarAz() const { return fStarAz*TMath::RadToDeg(); } + Double_t GetRawEl() const { return fRawEl*TMath::RadToDeg(); } + Double_t GetRawAz() const { return fRawAz*TMath::RadToDeg(); } + Double_t GetRawZd() const { return 90.-fRawEl*TMath::RadToDeg(); } ZdAz GetStarZdAz() const { return ZdAz(TMath::Pi()/2-fStarEl, fStarAz); } @@ -83,6 +119,6 @@ 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 AdjustEl(Double_t del) { fStarEl += del*TMath::DegToRad(); } + void AdjustAz(Double_t daz) { fStarAz += daz*TMath::DegToRad(); } void Adjust(const MBending &bend) @@ -105,20 +141,22 @@ 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; + TString str; + str.ReadLine(fin); + + Float_t v[4]; + sscanf(str.Data(), "%f %f %f %f", v, v+1, v+2, v+3); + + set.fStarAz = v[0]*TMath::DegToRad(); + set.fStarEl = v[1]*TMath::DegToRad(); + + set.fRawAz = v[2]*TMath::DegToRad(); + set.fRawEl = v[3]*TMath::DegToRad(); + + if (fin) + { + Double_t res, err; + res = set.GetResidual(&err); + cout << "Read: " << v[0] << " " << v[1] << " : " << v[2] << " " << v[3] << " : " << v[2]-v[0] << " " << v[3]-v[1] << " : " << res << " " << err << " " << err/res << endl; + } return fin; @@ -160,6 +198,7 @@ set.Adjust(bend); - Double_t err = 0.02; // [deg] = 1SE - Double_t res = set.GetResidual()/err; + Double_t err;// = 0.005; // [deg] = 0.25SE + Double_t res = set.GetResidual(&err); + res /= err; f += res*res; @@ -168,4 +207,5 @@ //f /= (fCoordinates.GetSize()-7)*(fCoordinates.GetSize()-7); //f /= fCoordinates.GetSize()*fCoordinates.GetSize(); + //cout << f << ": " << fCoordinates.GetSize() << endl; f /= fCoordinates.GetSize(); } @@ -187,24 +227,24 @@ TMarker mark0; - TMarker mark1; + //TMarker mark1; mark0.SetMarkerStyle(kStar); - mark1.SetMarkerStyle(kStar); - mark1.SetMarkerColor(kRed); + //mark1.SetMarkerStyle(kStar); + //mark1.SetMarkerColor(kRed); r0 /= 90; - r1 /= 90; - phi0 *= TMath::Pi()/180; - phi1 *= TMath::Pi()/180; - + //r1 /= 90; + phi0 *= TMath::DegToRad(); + //phi1 *= TMath::DegToRad(); + Double_t x0[3] = { r0*cos(phi0), r0*sin(phi0), 0}; - Double_t x1[3] = { r1*cos(phi1), r1*sin(phi1), 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); + //view->WCtoNDC(x1, y1); mark0.DrawMarker(y0[0], y0[1]); - mark1.DrawMarker(y1[0], y1[1]); + //mark1.DrawMarker(y1[0], y1[1]); } @@ -256,5 +296,5 @@ Double_t dp = p1-p0; - Double_t x0[3] = { r0*cos(p0*TMath::Pi()/180), r0*sin(p0*TMath::Pi()/180), 0}; + Double_t x0[3] = { r0*cos(p0*TMath::DegToRad()), r0*sin(p0*TMath::DegToRad()), 0}; for (double i=p0+10; iSetOptStat("emro"); + + TH1F hres1("Res1", " Residuals before correction ", fOriginal.GetSize()/3, 0, 0.3); + TH1F hres2("Res2", " Residuals after correction ", fOriginal.GetSize()/3, 0, 0.3); + TH1F hres3("Res3", " Residuals after backward correction ", fOriginal.GetSize()/3, 0, 0.3); hres1.SetXTitle("\\Delta [\\circ]"); @@ -505,6 +547,4 @@ cout << endl; - - TCanvas *c1; @@ -576,14 +616,20 @@ 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 << "before: " << Form("%.1f", hres1.GetMean()*60) << " \xb1 " << Form("%.1f", hres1.GetRMS()*60) << " arcmin" << endl; + cout << "after: " << Form("%.1f", hres2.GetMean()*60) << " \xb1 " << Form("%.1f", hres2.GetRMS()*60) << " arcmin" << endl; + cout << endl; + cout << "before: " << Form("%.1f", hres1.GetMean()*60*60/23.4) << " \xb1 " << Form("%.1f", hres1.GetRMS()*60*60/23.4) << " pix" << endl; + cout << "after: " << Form("%.1f", hres2.GetMean()*60*60/23.4) << " \xb1 " << Form("%.1f", hres2.GetRMS()*60*60/23.4) << " pix" << endl; + cout << "after: " << Form("%.1f", hres2.GetMean()*16384/360) << " \xb1 " << Form("%.1f", hres2.GetRMS()*16384/360) << " SE" << endl; + cout << endl; + cout << "backw: " << Form("%.1f", hres3.GetMean()*60) << " \xb1 " << Form("%.1f", hres3.GetRMS()*60) << " arcmin" << endl; cout << endl; // ± - + + before = hres1.GetMean()*16384/360; + after = hres2.GetMean()*16384/360; + backw = hres3.GetMean()*16384/360; + gStyle->SetOptStat(1110); @@ -615,5 +661,5 @@ 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; + cout << " Chi^2/NDF: " << f.GetChisquare() << "/" << f.GetNDF() << " = " << f.GetChisquare()/f.GetNDF() << endl; c1->cd(1); @@ -626,6 +672,6 @@ gPad->Modified(); gPad->Update(); - for (int i=0; icd(3); @@ -638,5 +684,5 @@ gPad->Update(); for (int i=0; i> set; // Read data from file [deg], it is stored in [rad] - if (!fin) break; @@ -687,13 +732,19 @@ { case kTbFit: - Fit(); - DisplayBending(); + { + Double_t before=0; + Double_t after=0; + Double_t backw=0; + Fit(before, after, backw); + DisplayBending(); + DisplayResult(before, after, backw); + } return kTRUE; case kTbLoad: - fBending.Load("bending_magic.txt"); + fBending.Load("bending_new.txt"); DisplayBending(); return kTRUE; case kTbSave: - fBending.Save("bending_magic.txt"); + fBending.Save("bending_new.txt"); return kTRUE; case kTbLoadStars: @@ -767,4 +818,16 @@ } + void DisplayResult(Double_t before, Double_t after, Double_t backw) + { + TGLabel *l1 = (TGLabel*)fLabel.At(3*MBending::GetNumPar()+1); + l1->SetText(Form("Before: %.1f +- %.1f SE", before)); + + TGLabel *l2 = (TGLabel*)fLabel.At(3*MBending::GetNumPar()+2); + l2->SetText(Form("After: %.1f +- %.1f SE", after)); + + TGLabel *l3 = (TGLabel*)fLabel.At(3*MBending::GetNumPar()+3); + l3->SetText(Form("Backw: %.1f +- %.1f SE", backw)); + } + public: ~MFit() @@ -866,4 +929,22 @@ fLabel.Add(l); + l = new TGLabel(grp1, ""); + l->SetTextJustify(kTextLeft); + grp1->AddFrame(l, hints5); + fList->Add(l); + fLabel.Add(l); + + l = new TGLabel(grp1, ""); + l->SetTextJustify(kTextLeft); + grp1->AddFrame(l, hints5); + fList->Add(l); + fLabel.Add(l); + + l = new TGLabel(grp1, ""); + l->SetTextJustify(kTextLeft); + grp1->AddFrame(l, hints5); + fList->Add(l); + fLabel.Add(l); + ((TGCheckButton*)fList->FindWidget(0))->SetState(kButtonDown); @@ -893,4 +974,5 @@ void gui() { + gErrorIgnoreLevel = kError; new MFit; }