Changeset 7217 for trunk/MagicSoft


Ignore:
Timestamp:
07/25/05 10:35:31 (19 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/Changelog

    r7216 r7217  
    2020
    2121                                                 -*-*- END OF LINE -*-*-
     22
     23 2005/07/25 Thomas Bretz
     24
     25   * mastro/MAstroSky2Local.[h,cc]:
     26     - allow setting a pointing deviation when calculating the rotation angle
     27
     28   * mhflux/MHDisp.[h,cc]:
     29     - added pot for a significance map
     30     - changed alignment range from 0.33/0.47 to 0.325/0.475
     31     - added fDeviation to calculation of Rotatio angle
     32     - removed some obsolete code
     33
     34   * mpointing/MPointingPos.[h,cc]:
     35     - added option to calculate Rotation angle using MPointingDev
     36
     37
     38
    2239 2005/07/22 Daniela Dorner
    2340

  • trunk/MagicSoft/Mars/mastro/MAstroSky2Local.cc

    r3568 r7217  
    129129// seen with an Alt/Az telescope.
    130130//
     131// dzd and daz is a pointing offset, which is subtracted from the
     132// calculated local coordinates correspoding to time, observatory
     133// and ra/dec.
     134//
    131135// For more information see MAstro::RotationAngle
    132136//
    133 Double_t MAstroSky2Local::RotationAngle(Double_t ra, Double_t dec) const
     137Double_t MAstroSky2Local::RotationAngle(Double_t ra, Double_t dec, Double_t dzd, Double_t daz) const
    134138{
    135139    TVector3 v;
    … …  
    137141    v *= *this;
    138142
    139     return MAstro::RotationAngle(ZZ(), XZ(), v.Theta(), v.Phi());
     143    return MAstro::RotationAngle(ZZ(), XZ(), v.Theta()-dzd, v.Phi()-daz);
    140144}

  • trunk/MagicSoft/Mars/mastro/MAstroSky2Local.h

    r3533 r7217  
    1818    MAstroSky2Local(const MTime &t, const MObservatory &obs);
    1919
    20     Double_t RotationAngle(Double_t ra, Double_t dec) const;
     20    Double_t RotationAngle(Double_t ra, Double_t dec, Double_t dzd=0, Double_t daz=0) const;
    2121
    2222    ClassDef(MAstroSky2Local, 1) // Rotation Matrix to convert sky coordinates to ideal local coordinates

  • trunk/MagicSoft/Mars/mhflux/MHDisp.cc

    r7208 r7217  
    5151#include "MParameters.h"
    5252
    53 #include "MHillasExt.h"
    54 #include "MHillasSrc.h"
     53#include "MHillas.h"
    5554#include "MSrcPosCam.h"
    5655#include "MPointingPos.h"
     56#include "MPointingDev.h"
    5757
    5858ClassImp(MHDisp);
    … …  
    6565//
    6666MHDisp::MHDisp(const char *name, const char *title)
    67     : fHilExt(0), fDisp(0), fSmearing(-1), fWobble(kFALSE)
     67    : fDisp(0), fDeviation(0), fSmearing(-1), fWobble(kFALSE)
    6868{
    6969    //
    … …  
    106106    }
    107107
    108     fHilExt = (MHillasExt*)plist->FindObject("MHillasExt");
    109     if (!fHilExt)
    110     {
    111         *fLog << err << "MHillasExt not found... aborting." << endl;
    112         return kFALSE;
    113     }
    114 
    115     fHilSrc = (MHillasSrc*)plist->FindObject("MHillasSrc");
    116     if (!fHilSrc)
    117     {
    118         *fLog << err << "MHillasSrc not found... aborting." << endl;
    119         return kFALSE;
    120     }
     108    fDeviation = (MPointingDev*)plist->FindObject("MPointingDev");
     109    if (!fDeviation)
     110        *fLog << warn << "MPointingDev not found... ignored." << endl;
    121111
    122112    MBinning binsx, binsy;
    … …  
    160150    Double_t rho = 0;
    161151    if (fTime && fObservatory && fPointPos)
    162         rho = fPointPos->RotationAngle(*fObservatory, *fTime);
     152        rho = fPointPos->RotationAngle(*fObservatory, *fTime, fDeviation);
    163153
    164154    // Vector of length disp in direction of shower
    … …  
    347337//
    348338// Return the mean signal in h around (0,0) in a distance between
    349 // 0.33 and 0.47deg
     339// 0.325 and 0.475deg
    350340//
    351341Double_t MHDisp::GetOffSignal(TH1 &h) const
    … …  
    360350        {
    361351            const Double_t d = TMath::Hypot(axex.GetBinCenter(x+1), axey.GetBinCenter(y+1));
    362             if (d>0.33 && d<0.47)
    363             //if (d>0.4 && d<0.8)
     352            if (d>0.325 && d<0.475)
    364353            {
    365354                sum += h.GetBinContent(x+1,y+1);
    … …  
    367356            }
    368357        }
    369 
    370358    return sum/cnt;
    371359}
    … …  
    415403            if (rmax<TMath::Hypot(px, py))
    416404                h2.SetBinContent(bin, 0);
     405            else
     406                if (h2.GetBinContent(bin)==0)
     407                    h2.SetBinContent(bin, 1e-10);
    417408        }
    418409}
    … …  
    427418    const TAxis &axey = *s.GetYaxis();
    428419
    429     for (int x=1; x<=axex.GetNbins(); x++)
    430         for (int y=1; y<=axey.GetNbins(); y++)
    431         {
    432             const Int_t  bin = s.GetBin(x,y);
    433 
    434             const Double_t sig = h1.GetBinContent(bin);
    435             const Double_t bg  = h2.GetBinContent(bin);
    436 
    437             const Double_t S = MMath::SignificanceLiMaSigned(sig, bg, TMath::Abs(scale));
    438 
     420    const Int_t n = TMath::Nint(0.2/axex.GetBinWidth(1));
     421
     422    for (int x=1+n; x<=axex.GetNbins()-n; x++)
     423        for (int y=1+n; y<=axey.GetNbins()-n; y++)
     424        {
     425            Double_t sig=0;
     426            Double_t bg =0;
     427
     428            for (int dx=-n; dx<=n; dx++)
     429                for (int dy=-n; dy<=n; dy++)
     430                {
     431                    if (TMath::Hypot((float)dx, (float)dy)>n)
     432                        continue;
     433
     434                    const Int_t  bin = s.GetBin(x+dx,y+dy);
     435
     436                    sig += h1.GetBinContent(bin);
     437                    bg  += h2.GetBinContent(bin);
     438                }
     439
     440            const Double_t S = sig>0 ? MMath::SignificanceLiMaSigned(sig, bg, TMath::Abs(scale)) : 0;
     441
     442            const Int_t bin = s.GetBin(x,y);
    439443            s.SetBinContent(bin, S);
    440444        }
    441445}
    442446
     447void MHDisp::Profile1D(const char *name, const TH2 &h) const
     448{
     449    if (!gPad)
     450        return;
     451
     452    TH1D *hrc = dynamic_cast<TH1D*>(gPad->FindObject(name));
     453    if (!hrc)
     454        return;
     455
     456    hrc->Reset();
     457
     458    const Double_t max = h.GetMaximum();
     459
     460    MBinning(50, -max*1.1, max*1.1).Apply(*hrc);
     461
     462    for (int x=1; x<=h.GetXaxis()->GetNbins(); x++)
     463        for (int y=1; y<=h.GetXaxis()->GetNbins(); y++)
     464        {
     465            const Int_t   bin  = h.GetBin(x,y);
     466            const Double_t sig = h.GetBinContent(bin);
     467            if (sig!=0)
     468                hrc->Fill(sig);
     469        }
     470
     471    gPad->SetLogy(hrc->GetMaximum()>0);
     472
     473    if (!fHistOff || hrc->GetRMS()<0.1)
     474        return;
     475
     476    // ---------- Fix mean ----------
     477    // TF1 *g = (TF1*)gROOT->GetFunction("gaus");
     478    // g->FixParameter(1, 0);
     479    // hrc->Fit("gaus", "BQI");
     480
     481    hrc->Fit("gaus", "QI");
     482
     483    TF1 *f = hrc->GetFunction("gaus");
     484    if (f)
     485    {
     486        f->SetLineWidth(1);
     487        f->SetLineColor(kBlue);
     488    }
     489}
    443490
    444491void MHDisp::Paint(Option_t *o)
    … …  
    460507    h1->SetContour(99);
    461508
     509    TH2 *hx=0;
     510
    462511    Double_t scale = 1;
    463512    if (fHistOff)
    … …  
    467516        TH2 *h=(TH2*)fHistOff->Project3D("yx_off");
    468517
    469         scale = -1;
    470 
    471518        const Double_t h1off = GetOffSignal(*h1);
    472519        const Double_t hoff  = GetOffSignal(fHistBg);
    … …  
    474521        scale = hoff==0 ? -1 : -h1off/hoff;
    475522
    476         const Bool_t sig = kFALSE;
    477 
    478         if (!sig)
    479             h1->Add(&fHistBg, scale);
    480         else
    481             MakeSignificance(*h1, *h1, fHistBg, scale);
     523        hx = (TH2*)pad->GetPad(4)->FindObject("Alpha_yx_sig");
     524        if (hx)
     525        {
     526            hx->SetContour(99);
     527            MakeSignificance(*hx, *h1, fHistBg, scale);
     528            MakeDot(*hx);
     529            MakeSymmetric(hx);
     530        }
     531
     532        h1->Add(&fHistBg, scale);
    482533
    483534        MakeDot(*h1);
     535        MakeSymmetric(h1);
    484536
    485537        delete h;
    486 
    487         MakeSymmetric(h1);
    488538    }
    489539
    … …  
    515565        func.SetParameter(0, h2->GetBinContent(1));
    516566        func.FixParameter(1, 0);
    517         func.SetParameter(2, 0.15);
     567        func.SetParameter(2, 0.12);
    518568        if (fHistOff)
    519569            func.FixParameter(3, 0);
    520         func.SetParameter(4, h2->GetBinContent(15));
     570        func.SetParameter(4, 0);//h2->GetBinContent(20));
    521571        h2->Fit(&func, "IQ", "", 0, 1.0);
    522572
    523573        h2->SetTitle(Form("P=(%.2f\\circ/%.2f\\circ) \\omega=%.2f\\circ f=%.2f", x0, y0, func.GetParameter(2), TMath::Abs(scale)));
    524574    }
     575
     576    pad->cd(5);
     577    if (h1)
     578        Profile1D("Exc", *h1);
     579
     580    pad->cd(6);
     581    if (hx)
     582        Profile1D("Sig", *hx);
    525583}
    526584
    … …  
    533591    AppendPad(o);
    534592
    535     TString name = Form("%s_1", pad->GetName());
    536     TPad *p = new TPad(name,name, 0.005, 0.005, 0.65, 0.995, col, 0, 0);
     593    // ----- Pad number 4 -----
     594    TString name = Form("%s_4", pad->GetName());
     595    TPad *p = new TPad(name,name, 0.5025, 0.3355, 0.995, 0.995, col, 0, 0);
     596    p->SetNumber(4);
     597    p->Draw();
     598    p->cd();
     599
     600    TH1 *h3 = fHist.Project3D("yx_sig");
     601    h3->SetTitle("Significance Map");
     602    h3->SetDirectory(NULL);
     603    h3->SetXTitle(fHist.GetXaxis()->GetTitle());
     604    h3->SetYTitle(fHist.GetYaxis()->GetTitle());
     605    h3->SetMinimum(0);
     606    h3->Draw("colz");
     607    h3->SetBit(kCanDelete);
     608
     609    if (fHistOff)
     610        GetCatalog()->Draw("white mirror same *");
     611
     612    // ----- Pad number 1 -----
     613    pad->cd();
     614    name = Form("%s_1", pad->GetName());
     615    p = new TPad(name,name, 0.005, 0.3355, 0.4975, 0.995, col, 0, 0);
    537616    p->SetNumber(1);
    538617    p->Draw();
    539618    p->cd();
    540619
    541     TH1 *h3 = fHist.Project3D("yx_on");
     620    h3 = fHist.Project3D("yx_on");
    542621    h3->SetTitle("Distribution of equivalent events");
    543622    h3->SetDirectory(NULL);
    … …  
    551630        GetCatalog()->Draw("white mirror same *");
    552631
     632    // ----- Pad number 2 -----
    553633    pad->cd();
    554634    name = Form("%s_2", pad->GetName());
    555     p = new TPad(name,name, 0.66, 0.005, 0.995, 0.5, col, 0, 0);
     635    p = new TPad(name,name, 0.005, 0.005, 0.2485, 0.3315, col, 0, 0);
    556636    p->SetNumber(2);
    557637    p->Draw();
    … …  
    559639    h3->Draw("surf3");
    560640
     641    // ----- Pad number 3 -----
    561642    pad->cd();
    562643    name = Form("%s_3", pad->GetName());
    563     p = new TPad(name,name, 0.66, 0.5, 0.995, 0.995, col, 0, 0);
     644    p = new TPad(name,name, 0.2525, 0.005, 0.4985, 0.3315, col, 0, 0);
    564645    p->SetNumber(3);
    565646    p->SetGrid();
    … …  
    577658    h->SetBit(kCanDelete);
    578659    h->Draw();
     660
     661    // ----- Pad number 5 -----
     662    pad->cd();
     663    name = Form("%s_5", pad->GetName());
     664    p = new TPad(name,name, 0.5025, 0.005, 0.7485, 0.3315, col, 0, 0);
     665    p->SetNumber(5);
     666    p->SetGrid();
     667    p->Draw();
     668    p->cd();
     669
     670    h3 = new TH1D("Exc", "Distribution of excess events", 10, -1, 1);
     671    h3->SetDirectory(NULL);
     672    h3->SetXTitle("N");
     673    h3->SetYTitle("Counts");
     674    h3->Draw();
     675    h3->SetBit(kCanDelete);
     676
     677    // ----- Pad number 6 -----
     678    pad->cd();
     679    name = Form("%s_6", pad->GetName());
     680    p = new TPad(name,name, 0.7525, 0.005, 0.9995, 0.3315, col, 0, 0);
     681    p->SetNumber(6);
     682    p->SetGrid();
     683    p->Draw();
     684    p->cd();
     685
     686    h3 = new TH1D("Sig", "Distribution of significances", 10, -1, 1);
     687    h3->SetDirectory(NULL);
     688    h3->SetXTitle("N");
     689    h3->SetYTitle("Counts");
     690    h3->Draw();
     691    h3->SetBit(kCanDelete);
    579692}
    580693

  • trunk/MagicSoft/Mars/mhflux/MHDisp.h

    r7193 r7217  
    1616class MHillasSrc;
    1717class MParameterD;
     18class MPointingDev;
    1819
    1920class MHDisp : public MHFalseSource
    2021{
    2122private:
    22     MHillasExt   *fHilExt;  //!
    23     MHillasSrc   *fHilSrc;  //!
    2423    MParameterD  *fDisp;    //!
     24    MPointingDev *fDeviation; //!
    2525
    2626    TH2D         fHistBg;
    27 
    2827    TH2D         fHistBg1;
    2928    TH2D         fHistBg2;
    … …  
    3938    Double_t DeltaPhiSrc(const TVector2 &v) const;
    4039
     40    void Profile1D(const char *name, const TH2 &h) const;
    4141    void Profile(TH1 &h2, const TH2 &h1, Axis_t x0=0, Axis_t y0=0) const;
    4242    void MakeSignificance(TH2 &s, const TH2 &h1, const TH2 &h2, const Double_t scale=1) const;

  • trunk/MagicSoft/Mars/mpointing/MPointingPos.cc

    r3666 r7217  
    4444#include "MTime.h"
    4545#include "MObservatory.h"
     46#include "MPointingDev.h"
    4647#include "MAstroSky2Local.h"
    4748
    … …  
    7172// For more information see MAstro::RotationAngle
    7273//
    73 Double_t MPointingPos::RotationAngle(const MObservatory &o, const MTime &t) const
     74Double_t MPointingPos::RotationAngle(const MObservatory &o, const MTime &t, const MPointingDev *dev) const
    7475{
    75     return MAstroSky2Local(t, o).RotationAngle(GetRaRad(), GetDecRad());
     76    return dev ?
     77        MAstroSky2Local(t, o).RotationAngle(GetRaRad(), GetDecRad(), dev->GetDevZdRad(), dev->GetDevAzRad()):
     78        MAstroSky2Local(t, o).RotationAngle(GetRaRad(), GetDecRad());
    7679}

  • trunk/MagicSoft/Mars/mpointing/MPointingPos.h

    r7205 r7217  
    1313class MTime;
    1414class MObservatory;
     15class MPointingDev;
    1516
    1617class MPointingPos : public MParContainer
    … …  
    5455
    5556    Double_t RotationAngle(const MObservatory &o) const;
    56     Double_t RotationAngle(const MObservatory &o, const MTime &t) const;
     57    Double_t RotationAngle(const MObservatory &o, const MTime &t, const MPointingDev *dev=0) const;
    5758    Double_t RotationAngle(const MObservatory &o, const MTime *t) const
    5859    {
Note: See TracChangeset for help on using the changeset viewer.