Ignore:
Timestamp:
03/11/03 13:52:48 (22 years ago)
Author:
tbretz
Message:
*** empty log message ***
File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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}
Note: See TracChangeset for help on using the changeset viewer.