Changeset 1810 for trunk


Ignore:
Timestamp:
03/11/03 13:52:48 (22 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Cosy
Files:
13 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Cosy/Changelog

    r1803 r1810  
    11                                                                  -*-*- 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
    247
    348 2003/03/02 - Daniela Dorner, Thomas Bretz (LaPalma):

  • trunk/MagicSoft/Cosy/base/MStar.h

    r1691 r1810  
    2121    void Set(Double_t mx, Double_t my) { fX=mx; fY=my; }
    2222
     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
    2339    ClassDef(MStar, 1)
    2440};

  • trunk/MagicSoft/Cosy/base/MStarList.cc

    r1691 r1810  
    11#include "MStarList.h"
     2
     3#include <iostream.h>
    24
    35void MStarList::RemoveTwins(Double_t radius)
    … …  
    1820            return;
    1921
     22        fStars.RemoveAt(idx);
     23
    2024        MStarListIter Next(this, *first, radius);
    21         Delete(idx);
    2225
    2326        MStar *pos;
    … …  
    2932            mx += pos->GetX();
    3033            my += pos->GetY();
    31             Delete(pos);
     34            fStars.Remove(pos);
    3235            cnt++;
    3336        }

  • trunk/MagicSoft/Cosy/base/MStarList.h

    r1760 r1810  
    3838
    3939    void RemoveTwins(Double_t radius);
     40
     41    void Sort() { fStars.Sort(); }
     42
     43    void Expand(int n) { fStars.Expand(n); }
    4044};
    4145

  • trunk/MagicSoft/Cosy/base/timer.cc

    r1793 r1810  
    2626//    fDiv = fmod((fMs+fSecs)/(60*60*24), 1.0);
    2727}
     28
     29void Timer::GetTimeval(struct timeval *tv) const
     30{
     31    tv->tv_sec  = fSecs;
     32    tv->tv_usec = fMs;
     33}
     34
    2835/*
    2936void Timer::Set(const long mjd)

  • trunk/MagicSoft/Cosy/base/timer.h

    r1760 r1810  
    3333    void SetTimer(const struct timeval *tv);
    3434
     35    void GetTimeval(struct timeval *tv) const;
     36
    3537    int GetSecs() { return fSecs; }
    3638

  • trunk/MagicSoft/Cosy/gui/MGAccuracy.cc

    r1760 r1810  
    195195    aaz *= d2r;
    196196
     197    const double el = TMath::Pi()/2-pzd;
     198
    197199    const double dphi2 = aaz/2.;
    198200    const double cos2  = cos(dphi2)*cos(dphi2);
    199201    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;
    201203
    202204    double dist = acos(d);
    … …  
    250252
    251253    UpdateCross(x, y);
    252     UpdateText(pos.Zd(), acc.Az(), acc.Zd());
     254    UpdateText(pos.Zd(), acc.Zd(), acc.Az());
    253255
    254256    SetModified();

  • trunk/MagicSoft/Cosy/gui/MGPngReader.cc

    r1802 r1810  
    33#include <fstream.h>    // ifstream
    44#include <iostream.h>   // cout
     5#include <iomanip.h>   // cout
    56
    67#include <TTimer.h>
    78
    89#include <TGMenu.h>
     10#include <TGLabel.h>
    911#include <TSystem.h>
    1012#include <TGSplitter.h>    // TGHorizontal3DLine
    … …  
    1315#include "MGImage.h"
    1416#include "MGCoordinates.h"
     17#include "MGDispCoordinates.h"
    1518
    1619#include "coord.h"
    … …  
    255258    fList->Add(fLimMag);
    256259
    257     fSao->SetLimitMag(8.0);
     260    fSao->SetLimitMag(7.0);
    258261
    259262    fInterpol = new MGPopupMenu(p);
    … …  
    314317    fList->Add(fCRaDec);
    315318
    316     fCZdAz = new MGCoordinates(this, kETypeZdAz);
     319    fCZdAz = new MGCoordinates(this, kETypeZdAz, kFALSE);
    317320    fCZdAz->Move(240+12+10, fMenu->GetDefaultHeight()+584);
    318321    AddFrame(fCZdAz);
    319322    fList->Add(fCZdAz);
    320323
     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
    321347    const Double_t pixsize = 23.4;
    322348
    … …  
    328354    fPixSize = new TGTextEntry(this, txt, IDM_kPixSize);
    329355    fPixSize->SetAlignment(kTextCenterX);
    330     fPixSize->Move(600, fMenu->GetDefaultHeight()+584);
     356    fPixSize->Move(547, fMenu->GetDefaultHeight()+617);
    331357    AddFrame(fPixSize);
    332358    fList->Add(fPixSize);
    … …  
    357383    MapSubwindows();
    358384    MapWindow();
     385
     386    //------------------------------------------------------------
     387    //    XY xy(3.819444, 24.05333);
     388    //    fCRaDec->SetCoordinates(xy);
     389    //    fRaDec->Set(xy.X()*360/24, xy.Y());
     390    //------------------------------------------------------------
    359391}
    360392
    361393MGPngReader::MGPngReader(MObservatory::LocationName_t obs)
    362 : TGMainFrame(gClient->GetRoot(), 768, 700), 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)
    363395{
    364396    fSao = new StarCatalog(obs);
    … …  
    811843}
    812844
    813 void MGPngReader::CalcTrackingError(MStarList &spots, MStarList &stars)
    814 {
     845ZdAz 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
     908void MGPngReader::CalcTrackingError(Leds &leds, MStarList &stars)
     909{
     910    const Int_t max = leds.GetEntries();
     911
    815912    if (stars.GetRealEntries() < 3)
    816913    {
    … …  
    819916    }
    820917
    821     if (spots.GetRealEntries() < 1)
     918    if (max < 1)
    822919    {
    823920        cout << "Sorry, less than 1 detected spot in FOV!" << endl;
    … …  
    825922    }
    826923
    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?
    830937
    831938    MStar *star;
    832     MStar *spot;
    833939    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()))
    856945        {
    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;
    869960        }
    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}
    905979
    906980
    … …  
    9451019
    9461020    MStarList spots;
     1021
     1022    /*
    9471023    if (fDisplay->IsEntryChecked(IDM_kStarguider))
    9481024        Filter2::Execute(spots, c);
    9491025    else
    950         if (fDisplay->IsEntryChecked(IDM_kFilter))
    951             Filter::Execute(c);
     1026     */
     1027    if (fDisplay->IsEntryChecked(IDM_kFilter))
     1028        Filter::Execute(c);
    9521029
    9531030    if (fDisplay->IsEntryChecked(IDM_kCaosFilter))
    … …  
    10311108    if (fDisplay->IsEntryChecked(IDM_kCatalog))
    10321109    {
     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
    10331135        DrawCircle(c, 0.5);
    10341136        DrawCircle(c, 1.0);
    10351137        DrawCircle(c, 1.5);
    10361138
    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());
    10471140
    10481141        fImage->DrawColImg(c, cimg);
    1049 
    1050         fCZdAz->SetCoordinates(fSao->GetZdAz());
    1051 
    1052         if (fDisplay->IsEntryChecked(IDM_kStarguider))
    1053             CalcTrackingError(spots, stars);
    10541142    }
    10551143    else
    … …  
    10671155    {
    10681156        const int dy = (int)sqrt(rpix*rpix-dx*dx);
    1069         img[cx+dx + (cy-dy)*768] = 0xd0;
    1070         img[cx+dx + (cy+dy)*768] = 0xd0;
    1071         img[cx-dy + (cy+dx)*768] = 0xd0;
    1072         img[cx+dy + (cy+dx)*768] = 0xd0;
     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;
    10731161    }
    10741162}

  • trunk/MagicSoft/Cosy/gui/MGPngReader.h

    r1802 r1810  
    1313#include "MGImage.h"
    1414
    15 #include <TH1.h>
    16 #include <TH2.h>
    17 #include <TGraph.h>
    18 #include <TCanvas.h>
     15#include "coord.h"
    1916
    20 class AltAz;
    21 class RaDec;
     17class TArrayF;
     18class TH1F;
     19class TH2F;
     20class TGraph;
    2221
    2322class TTimer;
    … …  
    8786    MGCoordinates *fCZdAz;
    8887
     88    MGCoordinates *fPZdAz;
     89
    8990    TGTextEntry   *fPixSize;
    9091
    … …  
    104105    void Toggle(MGPopupMenu *p, UInt_t id);
    105106    void GetCoordinates();
    106     void CalcTrackingError(MStarList &, MStarList &);
     107    void CalcTrackingError(Leds &, MStarList &);
     108    ZdAz TrackingError(TArrayF &alt, TArrayF &az, TArrayF &mag) const;
    107109
    108110    void InitHists();

  • trunk/MagicSoft/Cosy/main/MBending.cc

    r1808 r1810  
    107107        if (str=="ECEC")   fEcec = val;
    108108        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;
    111111
    112112        fin >> val;
    … …  
    143143    fout << "S   00   000000   000000  0000000" << endl;
    144144    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;
    161161    fout << "END" << endl;
    162162}
    … …  
    173173    p += CEC;
    174174
     175    //    const AltAz MAGIC(0, -fMagic1*tan(p.Alt())-fMagic2);
     176    //    p += MAGIC;
     177
    175178    const AltAz CRX(-fCrx*sin(p.Az()-p.Alt()),  fCrx*cos(p.Az()-p.Alt())/cos(p.Alt()));
    176179    const AltAz CRY(-fCry*cos(p.Az()-p.Alt()), -fCry*sin(p.Az()-p.Alt())/cos(p.Alt()));
    … …  
    183186    p += NRY;
    184187
    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;
    187190
    188191    const AltAz AW( fAw*sin(p.Az()), -fAw*cos(p.Az())*tan(p.Alt()));
    … …  
    191194    p += AN;
    192195
    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;
    196198
    197199    const AltAz CA(0, -fCa/cos(p.Alt()));
    … …  
    227229    p -= AW;
    228230
    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;
    231233
    232234    const AltAz NRY(fNry*cos(p.Alt()), -fNry*tan(p.Alt()));
    … …  
    285287        fEces =par[10]/kRad2Deg; // Elevation Centering Error (sin)
    286288    case 10:
    287         fCry  =par[9]/kRad2Deg; // Alt/Az Coude Displacement (E-W)
     289        fCry  =par[9] /kRad2Deg; // Alt/Az Coude Displacement (E-W)
    288290    case 9:
    289         fCrx  =par[8]/kRad2Deg; // Alt/Az Coude Displacement (N-S)
     291        fCrx  =par[8] /kRad2Deg; // Alt/Az Coude Displacement (N-S)
    290292    case 8:
    291         fNry  =par[7]/kRad2Deg; // Nasmyth rotator displacement, vertical
     293        fNry  =par[7] /kRad2Deg; // Nasmyth rotator displacement, vertical
    292294    case 7:
    293         fNrx  =par[6]/kRad2Deg; // Nasmyth rotator displacement, horizontan
     295        fNrx  =par[6] /kRad2Deg; // Nasmyth rotator displacement, horizontan
    294296    case 6:
    295         fAw   =par[5]/kRad2Deg; // Azimuth Axis Misalignment (E-W)
     297        fAw   =par[5] /kRad2Deg; // Azimuth Axis Misalignment (E-W)
    296298    case 5:
    297         fAn   =par[4]/kRad2Deg; // Azimuth Axis Misalignment (N-S)
     299        fAn   =par[4] /kRad2Deg; // Azimuth Axis Misalignment (N-S)
    298300    case 4:
    299         fCa   =par[3]/kRad2Deg; // Left-Right Collimation Error
     301        fCa   =par[3] /kRad2Deg; // Left-Right Collimation Error
    300302    case 3:
    301         fNpae =par[2]/kRad2Deg; // Az-El Nonperpendicularity
     303        fNpae =par[2] /kRad2Deg; // Az-El Nonperpendicularity
    302304    case 2:
    303         fIe   =par[1]/kRad2Deg; // Index Error in Elevation
     305        fIe   =par[1] /kRad2Deg; // Index Error in Elevation
    304306    case 1:
    305         fIa   =par[0]/kRad2Deg; // Index Error in Azimuth
     307        fIa   =par[0] /kRad2Deg; // Index Error in Azimuth
    306308    }
    307309}
    … …  
    312314    {
    313315    case 16:
    314         par[15]=fMagic2*kRad2Deg; //
     316        par[15]=kRad2Deg*fMagic2; //
    315317    case 15:
    316         par[14]=fMagic1*kRad2Deg; //
     318        par[14]=kRad2Deg*fMagic1; //
    317319    case 14:
    318         par[13]=fAcec*kRad2Deg; // Azimuth Centering Error (cos)
     320        par[13]=kRad2Deg*fAcec; // Azimuth Centering Error (cos)
    319321    case 13:
    320         par[12]=fEcec*kRad2Deg; // Elevation Centering Error (cos)
     322        par[12]=kRad2Deg*fEcec; // Elevation Centering Error (cos)
    321323    case 12:
    322         par[11]=fAces*kRad2Deg; // Azimuth Centering Error (sin)
     324        par[11]=kRad2Deg*fAces; // Azimuth Centering Error (sin)
    323325    case 11:
    324         par[10]=fEces*kRad2Deg; // Elevation Centering Error (sin)
     326        par[10]=kRad2Deg*fEces; // Elevation Centering Error (sin)
    325327    case 10:
    326         par[9]=fCry*kRad2Deg; // Alt/Az Coude Displacement (E-W)
     328        par[9] =kRad2Deg*fCry; // Alt/Az Coude Displacement (E-W)
    327329    case 9:
    328         par[8]=fCrx*kRad2Deg; // Alt/Az Coude Displacement (N-S)
     330        par[8] =kRad2Deg*fCrx; // Alt/Az Coude Displacement (N-S)
    329331    case 8:
    330         par[7]=fNry*kRad2Deg; // Nasmyth rotator displacement, vertical
     332        par[7] =kRad2Deg*fNry; // Nasmyth rotator displacement, vertical
    331333    case 7:
    332         par[6]=fNrx*kRad2Deg; // Nasmyth rotator displacement, horizontan
     334        par[6] =kRad2Deg*fNrx; // Nasmyth rotator displacement, horizontan
    333335    case 6:
    334         par[5]=fAw*kRad2Deg; // Azimuth Axis Misalignment (E-W)
     336        par[5] =kRad2Deg*fAw;  // Azimuth Axis Misalignment (E-W)
    335337    case 5:
    336         par[4]=fAn*kRad2Deg; // Azimuth Axis Misalignment (N-S)
     338        par[4] =kRad2Deg*fAn;  // Azimuth Axis Misalignment (N-S)
    337339    case 4:
    338         par[3]=fCa*kRad2Deg; // Left-Right Collimation Error
     340        par[3] =kRad2Deg*fCa;  // Left-Right Collimation Error
    339341    case 3:
    340         par[2]=fNpae*kRad2Deg; // Az-El Nonperpendicularity
     342        par[2] =kRad2Deg*fNpae; // Az-El Nonperpendicularity
    341343    case 2:
    342         par[1]=fIe*kRad2Deg; // Index Error in Elevation
     344        par[1] =kRad2Deg*fIe;  // Index Error in Elevation
    343345    case 1:
    344         par[0]=fIa*kRad2Deg; // Index Error in Azimuth
     346        par[0] =kRad2Deg*fIa;  // Index Error in Azimuth
    345347    }
    346348}
    … …  
    349351{
    350352    if (n<0)
    351         n = m.GetNumPars();
     353        n = 16;
    352354
    353355    Int_t ierflg = 0;
    … …  
    357359    case 16:
    358360        m.mnparm(15,"MAGIC2", fMagic2*kRad2Deg,  1, -360, 360, ierflg);
    359         // cout << "Init 3 CA:    " << fCa << endl;
    360361    case 15:
    361362        m.mnparm(14,"MAGIC1", fMagic1*kRad2Deg,  1, -360, 360, ierflg);
    362         // cout << "Init 3 CA:    " << fCa << endl;
    363363    case 14:
    364364        m.mnparm(13,"ACEC", fAcec*kRad2Deg,  1, -360, 360, ierflg);
    365         // cout << "Init 3 CA:    " << fCa << endl;
    366365    case 13:
    367366        m.mnparm(12,"ECEC", fEcec*kRad2Deg,  1, -360, 360, ierflg);
    368         // cout << "Init 3 CA:    " << fCa << endl;
    369367    case 12:
    370368        m.mnparm(11,"ACES", fAcec*kRad2Deg,  1, -360, 360, ierflg);
    371         // cout << "Init 3 CA:    " << fCa << endl;
    372369    case 11:
    373370        m.mnparm(10,"ECES", fEcec*kRad2Deg,  1, -360, 360, ierflg);
    374         // cout << "Init 3 CA:    " << fCa << endl;
    375371    case 10:
    376372        m.mnparm(9, "CRY",  fCry*kRad2Deg,  1, -360, 360, ierflg);
    377         // cout << "Init 3 CA:    " << fCa << endl;
    378373    case 9:
    379374        m.mnparm(8, "CRX",  fCrx*kRad2Deg,  1, -360, 360, ierflg);
    380         // cout << "Init 3 CA:    " << fCa << endl;
    381375    case 8:
    382376        m.mnparm(7, "NRY",  fNry*kRad2Deg,  1, -360, 360, ierflg);
    383         // cout << "Init 3 CA:    " << fCa << endl;
    384377    case 7:
    385378        m.mnparm(6, "NRX",  fNrx*kRad2Deg,  1, -360, 360, ierflg);
    386         // cout << "Init 3 CA:    " << fCa << endl;
    387379    case 6:
    388380        m.mnparm(5, "AW",   fAw*kRad2Deg,   1, -360, 360, ierflg);
    389         // cout << "Init 3 CA:    " << fCa << endl;
    390381    case 5:
    391382        m.mnparm(4, "AN",   fAn*kRad2Deg,   1, -360, 360, ierflg);
    392         // cout << "Init 3 CA:    " << fCa << endl;
    393383    case 4:
    394384        m.mnparm(3, "CA",   fCa*kRad2Deg,   1, -360, 360, ierflg);
    395         // cout << "Init 3 CA:    " << fCa << endl;
    396385    case 3:
    397386        m.mnparm(2, "NPAE", fNpae*kRad2Deg, 1, -360, 360, ierflg);
    398         // cout << "Init 2 NPAE:  " << fNpae << endl;
    399387    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);
    402389    case 1:
    403390        m.mnparm(0, "IA",   fIa*kRad2Deg,   1, -360, 360, ierflg);
    404         // cout << "Init 0 IA:    " << fIa << endl;
    405391    }
    406392}
    … …  
    408394void MBending::GetMinuitParameters(TMinuit &m, Int_t n=-1)
    409395{
    410     if (n<0)
     396    if (n<0 || n>m.GetNumPars())
    411397        n = m.GetNumPars();
    412398
    … …  
    465451    }
    466452}
    467 
     453/*
     454void 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*/
    468473void MBending::PrintMinuitParameters(TMinuit &m, Int_t n=-1) const
    469474{
    … …  
    471476        n = m.GetNumPars();
    472477
     478    cout << setprecision(3);
     479
    473480    Double_t par, err;
    474481
    … …  
    477484    case 16:
    478485        m.GetParameter(15, par, err);
    479         cout << " 15 MAGIC2: " << par << " +- " << err << endl;
     486        cout << " 15 MAGIC2: " << setw(8) << par << " +- " << setw(4) << err << endl;
    480487    case 15:
    481488        m.GetParameter(14, par, err);
    482         cout << " 14 MAGIC1: " << par << " +- " << err << endl;
     489        cout << " 14 MAGIC1: " << setw(8) << par << " +- " << setw(4) << err << endl;
    483490    case 14:
    484491        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;
    486493    case 13:
    487494        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;
    489496    case 12:
    490497        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;
    492499    case 11:
    493500        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;
    495502    case 10:
    496503        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;
    498505    case 9:
    499506        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;
    501508    case 8:
    502509        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;
    504511    case 7:
    505512        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;
    507514    case 6:
    508515        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;
    510517    case 5:
    511518        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;
    513520    case 4:
    514521        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;
    516523    case 3:
    517524        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;
    519526    case 2:
    520527        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;
    522529    case 1:
    523530        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  
    11#include <fstream.h>
    22#include <iostream.h>
     3#include <iomanip.h>
    34
    45#include <TF1.h>
    … …  
    2122// Sekans = 1/cos
    2223//
    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 }
    3124
    3225class Set : public TObject
    … …  
    5447    Double_t GetResidual() const
    5548    {
    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;
    6151
    6252        Double_t dphi2 = daz/2.;
    6353        Double_t cos2  = cos(dphi2)*cos(dphi2);
    6454        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);
    6858
    6959        return dist * 180 / TMath::Pi();
    … …  
    7969
    8070    Double_t GetDEl() const     { return (fRawEl-fStarEl)*180/TMath::Pi(); }
     71    Double_t GetDZd() const     { return -GetDEl(); }
    8172    Double_t GetDAz() const     { return (fRawAz-fStarAz)*180/TMath::Pi(); }
    8273    Double_t GetStarEl() const  { return fStarEl*180/TMath::Pi(); }
    … …  
    10293        fStarAz = p.Az();
    10394    }
    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();
    112100    }
    113101};
    … …  
    116104{
    117105    Double_t v[4];
     106    fin >> v[0];
     107    fin >> v[1];
    118108    fin >> v[2];
    119109    fin >> v[3];
    120     fin >> v[0];
    121     fin >> v[1];
    122110
    123111    set.fStarAz = v[0]*TMath::Pi()/180;
    … …  
    140128
    141129    MBending bend;
    142     bend.SetParameters(par/*, npar*/); // Set Parameters [deg] to MBending
     130    bend.SetParameters(par); // Set Parameters [deg] to MBending
    143131
    144132    for (int i=0; i<gCoordinates.GetSize(); i++)
    … …  
    146134        Set set = *(Set*)gCoordinates.At(i);
    147135
    148         //set.Adjust(bend);
    149         set.Adjust(bend, MyAdjust);
     136        set.Adjust(bend);
    150137
    151138        Double_t err = 1;
    … …  
    200187        return;
    201188    }
    202 
     189    /*
    203190    if (r0<0)
    204191    {
    … …  
    212199    }
    213200
    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;
    217206    */
    218 
    219207    TLine line;
    220208    line.SetLineWidth(2);
    … …  
    268256    Double_t phi1 = set.GetStarAz();
    269257
     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
    270278    if (scale>0)
    271279    {
    272280        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;
    276283        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;
    280286    }
    281287
    … …  
    290296    gCoordinates.SetOwner();
    291297
    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 
    295298    TH1F hres1("Res1", "  Residuals before correction     ", 10, 0, 180);
    296299    TH1F hres2("Res2", "  Residuals after correction     ",  10, 0,   1);
    297300
    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 
    307301    hres1.SetXTitle("\\Delta [\\circ]");
    308302    hres1.SetYTitle("Counts");
    … …  
    313307    TGraph gdaz;
    314308    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");
    317319
    318320    while (fin && fin.get()!='\n');
    … …  
    346348    bending.SetMinuitParameters(minuit, 16); // Init Parameters [deg]
    347349
    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
    356359    minuit.FixParameter(8);  // CRX
    357360    minuit.FixParameter(9);  // CRY
    358    // minuit.FixParameter(10); // ECES
    359    //  minuit.FixParameter(11); // ACES
    360    // minuit.FixParameter(12); // ECEC
    361    //  minuit.FixParameter(13); // ACEC
    362     // minuit.FixParameter(14); // MAGIC1
     361    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
    363366    minuit.FixParameter(15); // MAGIC2
    364367
     368
    365369    //minuit.Command("SHOW PARAMETERS");
     370    //minuit.Command("SHOW LIMITS");
     371    cout << endl;
     372
     373    Int_t ierflg = 0;
    366374    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    //
    381383    bending.GetMinuitParameters(minuit);
    382384    bending.PrintMinuitParameters(minuit);
    383385    bending.Save("bending_magic.txt");
    384386
     387
     388    //
     389    // Make a copy of all list entries
     390    //
    385391    TList list;
    386392    list.SetOwner();
    387 
    388393    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++)
    389400    {
    390401        Set &set0 = *(Set*)gCoordinates.At(i);
    391402
    392         list.Add(new Set(set0));
     403        ZdAz za = set0.GetStarZdAz()*kRad2Deg;
    393404
    394405        hres1.Fill(set0.GetResidual());
    395         //set0.Adjust(bending);
    396         set0.Adjust(bending, MyAdjust);
     406        set0.Adjust(bending);
    397407        hres2.Fill(set0.GetResidual());
    398408
    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;
    401412
    402413        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());
    405418        j++;
    406419    }
    407420
     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    //
    408433    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;
    409462    cout << "Spead before: " << hres1.GetMean() << "deg  +-  " << hres1.GetRMS() << "deg" << endl;
    410463    cout << "Spead after:  " << hres2.GetMean() << "deg  +-  " << hres2.GetRMS() << "deg" << endl;
    … …  
    414467    cout << endl;
    415468
     469
     470
    416471    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*");
    418476
    419477    c1->cd(1);
    420     hdaz.DrawCopy("cont");
    421     gdaz.DrawClone("*");
     478    gaz.DrawClone("A*");
     479
     480    c1->cd(3);
     481    gdzd.DrawClone("A*");
    422482
    423483    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);
    426490
    427491    c1->cd(2);
    … …  
    429493    hres1.DrawCopy();
    430494
    431     c1->cd(5);
     495    c1->cd(4);
    432496    hres2.SetLineColor(kBlue);
    433497    hres2.DrawCopy();
    434498
    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());
    452518
    453519}

  • trunk/MagicSoft/Cosy/videodev/CaosFilter.cc

    r1802 r1810  
    345345    }
    346346
     347    RemoveTwins(leds, 5);
     348}
     349
     350void CaosFilter::RemoveTwins(Leds &leds, Double_t radius)
     351{
    347352    TIter Next1(&leds);
    348353
    … …  
    368373            const Double_t dy = y2-y1;
    369374
    370             if (dx*dx+dy*dy<5*5)
     375            if (dx*dx+dy*dy<radius*radius)
    371376            {
    372377                // FIXME: Interpolation

  • trunk/MagicSoft/Cosy/videodev/CaosFilter.h

    r1803 r1810  
    4949    static void Execute(byte *img, Leds &leds, Double_t conv);
    5050
    51     static void  FilterLeds(Leds &leds);
     51    static void FilterLeds(Leds &leds);
     52    static void RemoveTwins(Leds &leds, Double_t radius);
    5253
    5354    ClassDef(CaosFilter, 0)
Note: See TracChangeset for help on using the changeset viewer.